1.

Introduction

A Sun-like star is much brighter (typically ${10}^8$ to ${10}^{10}$ times) than an Earth-like planet in its habitable zone.¹ Moreover, at a distance of 10 pc, the star and planets in its habitable zone are separated by around 0.1 arc sec. Thus it is difficult to separate the planet light from that of the star in the image. There are two main solutions to the challenge of imaging objects in close proximity to much brighter ones. First, one can use a coronagraph,² which is a device inside the telescope to block the starlight from reaching the image plane. Second, one can use a starshade,^3,4 which is a large screen flying on a separate spacecraft positioned between the telescope and the star being observed to suppress the starlight before it enters the telescope. In many ways, coronagraphs and starshades are complementary. Coronagraphs are efficient for high-contrast surveys, because it is easy to point the instrument to different targets. However, they result in a lower optical throughput relative to a starshade, have the difficulty for designing a coronagraph due to “exotic pupils,” and are very sensitive to wavefront perturbations. Even small aberrations introduce bright speckles, which influence the instrument’s ability for exoplanet observations.⁵ Thus it imposes many challenging requirements on the telescope and instruments to design coronagraphs with $\sim {10}^{-10}$ starlight suppression; two mission concepts under study are Habitable Exoplanet Observatory (HabEx)⁶ and Large UV/Optical/Infrared Surveyor.⁷ In comparison, starshades are good at deep imaging and spectroscopic characterization. They are not sensitive to wavefront errors and can be designed to operate over a large bandpass. The total throughput is high since the starshade does not require any internal masking of the optical beam, which makes a starshade an excellent option for deep spectroscopy, especially at small inner working angles (IWA). However, one disadvantage of a starshade is the time it takes to slew the starshade to realign it with different target stars. The starshade’s ability to efficiently suppress the on-axis starlight while maintaining high throughput makes it an excellent tool for exploring the habitability of exoplanets. A recently studied potential mission is the Starshade Rendezvous Mission: a starshade that will work with the Nancy Grace Roman Space Telescope (previously called WFIRST).⁸ In the mission, the starshade is launched separately and rendezvous with the telescope in orbit. Starshades are also baselined for the HabEx mission concept.⁶

Starshades are a new technology, still in development. Coronagraphs, however, have been used on ground-based telescopes for decades; even the Hubble Space Telescope has a rudimentary coronagraph. Available research on high-contrast imaging is mostly about image processing and signal detection for coronagraph observations, which focus on alleviating the influence of quasi-static speckles. However, quasi-static speckles are not present in starshade images, so the emphasis is on the starshade’s error sources. The dominant sources of noise in starshade images are sunlight scattering off the starshade edges⁹ and unsuppressed starlight caused by errors in the starshade shape.¹⁰ The scattered sunlight is confined to two extended lobes perpendicular to the direction of the Sun and is constant during observations.¹¹ Its stability means that it could potentially be calibrated out, just adding photon noise to a small region around the starshade. Manufacturing and deployment errors and thermal deformations can distort the starshade shape and will produce bright spots in images that are hard to distinguish from a real planet signal. One example shape error we examine in this study is the truncation of the tips of starshade petals. Additional sources of noise in starshade images come from misalignment of the starshade and telescope, detector noise, and exozodiacal dust.¹²

Due to the difference of the noise properties in coronagraph images and starshade images, previous work on coronagraphs loses its utility on starshade images, which motivates the investigation of new techniques for starshade images. In this paper, we focus on the impact caused by errors in the starshade shape and exozodiacal dust and present an automatic detection algorithm, the generalized likelihood ratio test (GLRT), to provide robust detection on low-signal images in the presence of shape errors. We have described the GLRT model and its preliminary results for simulated images with starshade shape errors, dynamics, and detector noise in Ref. 13. We will review this detection method and introduce an iterative process to detect a planet in the presence of significant exozodiacal dust. This work focuses on signal detection without the need for postprocessing [e.g., point spread function (PSF) subtraction]. Post-processing may improve the detection performance but could also complicate the data analysis process and risks introducing artifacts into or removing part of the planetary signal. We believe demonstrating our signal detection method on unprocessed images strengthens the argument for the efficiency of our method. Detailed investigation on post-processing is beyond the scope of the paper.

In Sec. 2, we describe the image simulation process used to generate the test set for our planet detection methods. Section 3 presents the GLRT detection method and represents the bulk of this work. An iterative approach of GLRT is presented in Sec. 4 to tackle the problem of exozodiacal dust. We end by summarizing our work.

2.

Image Generation

We briefly summarize the image generation process outlined in Fig. 1. A more detailed treatment can be found in Ref. 14. The input for the image generation process is a model of the solar system viewed face-on from 10 pc away developed as part of the Haystacks Project.¹⁵ This model contains multi-wavelength image slices, covering the range from 0.3 to $2.5\mu \mathrm{m}$. Each image centers on the star, extending to 50 AU from it. The star and planets are represented as single-pixel sources and the model also includes interplanetary dust and astrophysical background sources. The pixel values in the image slices are spectral flux densities.

Fig. 1

Diagram of the image simulation process with starshade system illustration (not drawn to scale). The illustration, i.e., light sources, a starshade and a telescope, on the bottom is aligned with the description above. ROI is the area defined at the simulation input, i.e., the astronomical scene, beyond which the incoming light is considered a plane wave and is not diffracted by the starshade.

The optical model to calculate the starshade diffraction uses Fresnel diffraction theory¹⁶ to propagate light past the starshade. However, calculating the propagation of each pixel separately in Haystacks is computationally expensive. Starshades have a noticeable influence on light propagation only for a small area in close proximity to the starshade, which we call the region of influence (ROI). The ROI is defined at the input plane of the simulation and is the angular separation of a source, beyond which the incoming light is considered a plane wave and is not diffracted by the starshade. The result for an off-axis light source outside the ROI is close to the result for the same light source as if there were no starshade. Thus we only use Fresnel propagation inside the ROI and simply convolve point sources outside the ROI with the telescope’s PSF. The starshade we use in this paper is 13 m in radius and has 16 petals, shown in Fig. 2.

Fig. 2

The starshade we use in this paper, which is 13 m in radius and has 16 petals and designed for the Starshade Rendezvous Mission.⁸

The image simulation incorporates the main factors that influence the image of a realistic system: a matrix represents the starshade outline; defects are added to the shape by adjusting this matrix; the alignment between the starshade and telescope is accounted for by adding time-dependent formation flying dynamics to the diffraction calculation; and a detector model for the Roman Telescope is used, which includes clock-induced charge, dark current, degradation during lifetime, polarization losses, and read noise.¹⁷

Planets are extremely faint, so their signals can be weak compared to the readout noise. To tackle this problem, the Roman Telescope uses an electron-multiplication charge-coupled device (EMCCD),¹⁷ which amplifies the signal in an electron-multiplication register. This process reduces the effective readout noise to less than one electron.¹⁸ EMCCDs introduce an additional noise source: the multiplicative noise associated with the amplification process, which can be eliminated using the detector in photon counting (PC) mode. PC mode applies a chosen threshold to the number of electrons generated in each pixel and yields a value of one if the number of electrons is larger than the threshold, zero otherwise. It is a binary process that does not measure the exact number of photons but rather the presence of photons. As PC mode cannot distinguish the event of one photon from the event of more than one photon, the exposure time is short enough so that the expected number of photons in any pixel of the detector is $<1$. In this way, photons are not wasted. In this work, we choose the integration time for each PC image to be 1 s so that the maximum photon rate on the detector (looking at an Earth-like planet from 10 pc) is around 0.1 photon per second per pixel.

Typically, the “binary” PC images are not used directly but rather are added together to create a final image, which we will call a co-added image. In this way, the co-added images have a large enough dynamic range so the different intensities of the sources can be well-reflected. Most of the images shown in this paper are co-added final images, unless otherwise stated. The number of PC images to combine into one co-added image is denoted by $N_{\mathrm{im}}$. In this work, we use $N_{\mathrm{im}}=2000$, if not specified, otherwise; guidance on how to choose $N_{\mathrm{im}}$ is provided in Sec. 3.6. The imaging field-of-view diameter is 9 arc sec and each pixel is $0.021\mathrm{arc}\sec \times 0.021\mathrm{arc}\sec$. To visualize the detector’s performance, we use Monte Carlo to calculate the probability density functions (PDFs) of the photon counts for different ground-truth photon fluxes in a co-added image from 2000 PC images, shown in Fig. 3(a). Four of the PDF’s are plotted in Fig. 3(b); the parameters in the simulations are listed in Table 1. A photon rate of $0.1\text{photon}/\mathrm{s}$ (our expected maximum rate) is well within the linear response regime.

Fig. 3

PDFs of the photon counts for different ground-truth photon fluxes. (a). Each row is the PDF of photon counts observed in one pixel in a co-added image from 2000 PC images corresponding to the photon flux in that pixel. (b) Example PDF’s of the photon counts for four different ground-truth photon fluxes, which correspond to the ones in (a). The PDF of photon counts observed in one pixel in a co-added image from 2000 PC images corresponding to the photon flux 0.01, 0.1, 1, and $5\text{count}/\mathrm{s}$ in the pixel.

Table 1

Parameters for simulation of solar system as viewed from 10 pc.

Figure 4 shows one wavelength slice of a Haystacks model, along with a few results from different stages in our simulation. Figure 4(a) shows the input of our simulation, a discretized astrophysical scene of our solar system as viewed from 10 pc. The pixel size is small enough to make sure that the Haystacks models are high-fidelity spatially. The Haystacks scene includes the star, planets, interplanetary dust, and astrophysical background sources. As described in Fig. 1, the light sources, i.e., Fig. 4(a), are diffracted at the starshade plane and then propagated to the telescope’s pupil plane. A simple (aberration-free) telescope model propagates the pupil plane to the detector plane, the result of which is shown in Fig. 4(b). Then Fig. 4(b) is processed with the detector model to generate one single PC image, which is shown in Fig. 4(c). By combing 2000 PC images, we get the co-added image [Fig. 4(d)]. The dark region in the center is the starlight suppression effect of the starshade. The brighter ring is the exozodiacal dust outside the starshade’s IWA. Mars is too dim to be seen and Earth is of similar brightness as the exozodiacal dust. Venus is bright enough to be seen. In the following sections, we propose a detection and characterization method of planetary signals for co-added images at $\lambda =0.55\mu \mathrm{m}$ and with a bandwidth of $\lambda =0.12\mu \mathrm{m}$.

Fig. 4

The solar system as viewed from 10 pc—and its results from different stages in our simulation. They are zoomed-in views of the images. (a) The original astronomical scene from Haystacks project¹⁵ at $\lambda =0.55\mu \mathrm{m}$ in log10 scale flux (the central star is made ${10}^{10}$ times dimmer to reveal fainter features). (b) A simulated result before including the detector noise, with a perfect starshade at $\lambda =0.55\mu \mathrm{m}$ and with a bandwidth of $\lambda =0.12\mu \mathrm{m}$. (c) One PC image. (d) the co-added image from 2000 PC images.

3.

Generalized Likelihood Ratio Test

This section presents the GLRT as an automatic detection method for starshade missions. Our work is motivated by the lack of any previous investigation into signal detection in starshade images. We begin the section by reviewing past work on signal detection for direct imaging, most of which are specialized to coronagraphic images.

The biggest challenge for coronagraphic images is the noise floor set by quasi-static speckles. Different observing techniques and post-processing methods have been developed to try to attenuate the speckles before attempting detection. They are all based on differential imaging, which consists of estimating the star-only coronagraphic image and subtracting it from the science images (also called speckle subtraction technique). Differential imaging relies on specific observation strategies such as angular diversity (ADI),¹⁹ spectral diversity,²⁰ or multiple reference star images²¹ to generate the differential signal. The various speckle subtraction methods developed for coronagraphic images may serve to improve detection in starshade images, but it is beyond the scope of this work to include them. We focus on planet detection in images that have not been post-processed and leave that work to be done in the future.

In our work, we do not post-process the images (as in speckle subtraction), so we now move to detection methods. Braems and Kasdin^22,23 used hypothesis testing for the detection. The unknown parameters such as the planets’ positions and intensities can be removed by marginalizing the probability²² or using worst-case values.²³ However, the choice of priors or using the worst-case values is an open question. Kasdin and Braems²³ assumed a known constant background; in our method, we will estimate the background with a maximum likelihood. Mawet et al.²⁴ (SNRt map) essentially tested if the intensity of a single test resolution is different from the other resolution elements in a $1\lambda /D$-wide annulus at the same radius. They used small sample statistics to address the problem of the statistical significance of detection when few realizations of spatial speckles versus azimuth are present. Cantalloube et al.²⁵ (ANDROMEDA) made a signal template considering over/self subtraction caused by ADI rather than directly using a theoretical PSF template when calculating maximum likelihood estimation (MLE). Ruffio et al.²⁶ (FMMF) also used Signal-to-noise ratio (SNR)-based methods. They include Karhunen-Lòeve Image Projection (KLIP)-induced distortion in the match-filter template to estimate the signal intensity. The standard deviation is calculated at each pixel while masking a disk with a 5-pixel radius centered on that pixel from the annulus to prevent a planet biasing its own SNR. Flasseur et al.²⁷ (PACO) used GLRT and used the full-covariance rather than assume independent and identically distributed in the time dimension for the different frames. However, they assumed that the covariance is the same under $H0$, $H1$ and thus calculate only one covariance. Moreover, this covariance is not derived from the Gaussian equation together with the signal intensity, and thus the result is not guaranteed to be MLE. Pairet et al.²⁸ (STIM map) proposed STIM map, which used a modified Rician distribution, which is a heavy tail distribution compared to Gaussian, to model the speckles. Dahlqvist et al.²⁹ (RSM map) used a Markov process to model the same pixel throughout the different images. They claimed the residual quasi-static speckles after ADI can be characterized by their mean and variance. Recent studies also took the detection problem as a binary classification problem, using random-forest classifiers (SODIRF) and deep neural networks (SODINN).³⁰ However, these machine learning methods need very large training sets, which are difficult to generate for unknown planet signals, and require a heuristic tuning of hyper-parameters.

In this study, we also take the detection problem as a “hypothesis testing” problem and use a GLRT. The null hypothesis ($H_0$) is that there is no planet; the alternative hypothesis ($H_1$) is that there is a planet. We compare the posterior probabilities under the two hypotheses to decide whether to reject the null. Instead of marginalization, we use MLE to first estimate the unknown parameters (intensity, position, and standard deviation of the noise) and then use them to calculate the likelihood ratio. Using this ratio, even when the maximum likelihood under the alternative hypothesis is low, we may still have a strong detection as long as the ratio is high (i.e., the pattern is much less likely from pure noise). We first introduced the GLRT model and its preliminary results in Ref. 13. In this section, we briefly describe the method, show recent improvements, and present results.

3.4.

Detection in the Whole Image

We have built a model and corresponding detection and estimation method for a search area that is the size of the PSF core. However, the potential planets’ locations are unknown, so to perform detection in the whole image, we will traverse the whole image using the method outlined in Sec. 3.3. For each pixel, we use the detection area centered at this pixel, that is to say, we test $H_1$: there is a planet centered at this pixel, against $H_0$: there is only constant background there. After calculating for all the pixels in the image, we choose an FA rate and apply its corresponding threshold. Examples are illustrated in Fig. 5. As the PSF changes with the distance from the starshade, our detection algorithm uses a library of reference PSFs.

Fig. 5

Example of the GLRT detection on an image with perfect starshade and without exozodiacal dust. (a) Noiseless image: it should be taken as ground truth for image processing. (b) Image with detector noise: examples of search areas are also shown. The white box is centered on pixel (7,12), which is marked by the white asterisk. The pixel values in the white box form data ${\mathit{x}}_{7,12}$ in Eq. (3). The magenta box forms data ${\mathit{x}}_{8,14}$ and is the case where the search area is at the edge of the PSF. (c) The $T$ values from Eq. (24) in each pixel. (d) FA rate map. (e) binary detection image after thresholding. We apply a threshold of 0.7354 to the $T$ map, which results in 0.4 FA rate. The white circles are the minimal bounding circles used to estimate the planet position. (f) The relationship between threshold and FA rate from Eq. (25).

When the pixel is at the boundary of the PSF in the image, [one example is shown as a magenta asterisk in Fig. 5(b)], the detection area only contains part of the planet which is not centered at the pixel, so neither $H_1$ nor $H_0$ is true and the MLE of planet intensity can be negative. Thus we set those negative estimates as zero and thus $T({\mathit{x}}_{i,j})=0$.

After thresholding, we get a binary image. Generally speaking, some pixels next to the signal center will also be detected. Thus to estimate the position of the planet, we first find the convex hulls in the thresholded image and find the minimal circle bounding each convex hull. The center of the circles will be the estimates for the planets’ positions. One example is shown in Fig. 5(e). The estimated planet intensity is the MLE of $I$ at the estimated planet position. As we have shown in our previous work,¹⁴ the PSF changes with the distance away from the starshade center. If we use only one PSF template in GLRT, normally the one without a starshade, we can have a higher FA rate, worse position estimation, and worse intensity estimation. Thus we need to have a library of PSFs at difference distances from the starshade center for the GLRT model. In this paper, we define $\text{intensity error}=\frac{(\text{estimated intensity}-\text{real intensity})}{\text{real intensity}}$ and report its value for all cases.

3.5.

Results

Two detection and estimation examples are presented in Figs. 5 and 6. In Fig. 5, the detection process is demonstrated step by step with an image with perfect starshade. In Fig. 6, we demonstrate the GLRT’s ability to distinguish real signals from fake ones by providing the detection result for an image with a starshade with truncated tips. Compared to the perfect starshade, the tip of one of this starshade’s petals is truncated by 6.5 mm and resulting in a tip width of 48.3035 mm (the radius of the designed starshade is 13 m and the petal tip width is 48.3216 mm). We refer to this starshade as the clipped starshade. This defect causes a bright spot in the image, which could be mistaken for a planet. We still use the PSF templates of the perfect starshade for this clipped starshade case. GLRT successfully detects Venus and Earth from the fake signal. The errors of intensity and position estimation for these two examples are presented in Table. 2. The fake signal is close to Venus, so the intensity estimation of Venus is degraded for the clipped starshade case. The pixel size of the images is 0.021 arc sec, so the position estimation is accurate to the pixel level.

Fig. 6

Example of the GLRT detection on an image with clipped starshade. (a) Noiseless image: the defect on the starshade causes a bright spot that resembles a planet; (b) image with detector noise applied; (c) the $T$ values from Eq. (24) in each pixel; (d) FA rate map; and (e) binary detection image after thresholding. We apply a threshold of 0.7354 to the $T$ map, which results in 0.4 FA rate. GLRT successfully detects Venus and Earth and ignores the fake signal in the image.

Table 2

Intensity estimation error and position estimation error comparison between results in images using different starshades.

Starshade	Intensity error for Venus (%)	Intensity error for Earth (%)	Position error for Venus (arc sec)	Position error for Earth (arc sec)
Perfect starshade	1.5	1.2	$3\times {10}^{-3}$	$9.5\times {10}^{-3}$
Clipped starshade	20.5	4.1	$3\times {10}^{-3}$	$9.5\times {10}^{-}3$

We also calculate the receiver operating characteristic (ROC) curves for Venus and Earth shown in Fig. 7, where we compare the performance of GLRT for co-added images with different numbers of total images, which is denoted as $N_{\mathrm{im}}$. Equations (25) and (26) give the theoretical FA rate and true positive (TP) rate under a Gaussian assumption and are therefore only an approximation. Moreover, as shown in Eq. (27), the calculation of the TP rate needs the value of the true variance, which we need to estimate. Thus to more accurately demonstrate the detection’s performance, we use a Monte Carlo simulation. To calculate the ROC curves, we apply GLRT to get the FA rate map for a set of different thresholds and record if it results in a detection or missed detection of Earth, Venus, and a background pixel. We run a large number of trials, which is denoted as $n_{\text{trials}}$ and record the ratio of detection of Earth and Venus as the TP rate for Earth and Venus and record the ratio of detection of the background pixel as the false positive rate. The CI of a proportion $\widehat{p}$ is³⁵

Eq. (28)

$$(\widehat{p}-z\sqrt{\frac{\widehat{p}(1-\widehat{p})}{n_{\text{trials}}}},\widehat{p}+z\sqrt{\frac{\widehat{p}(1-\widehat{p})}{n_{\text{trials}}}}).$$

Thus for a point $({\widehat{p}}_{\text{false}},{\widehat{p}}_{\text{positive}})$ on a ROC curve, we take the CI using the two points:

Eq. (29)

$$({\widehat{p}}_{\text{false}}-z\sqrt{\frac{{\widehat{p}}_{\text{false}}(1-{\widehat{p}}_{\text{false}})}{n_{\text{trials}}}},{\widehat{p}}_{\text{positive}}+z\sqrt{\frac{{\widehat{p}}_{\text{positive}}(1-{\widehat{p}}_{\text{positive}})}{n_{\text{trials}}}}),$$

and

Eq. (30)

$$({\widehat{p}}_{\text{false}}+z\sqrt{\frac{{\widehat{p}}_{\text{false}}(1-{\widehat{p}}_{\text{false}})}{n_{\text{trials}}}},{\widehat{p}}_{\text{positive}}-z\sqrt{\frac{{\widehat{p}}_{\text{positive}}(1-{\widehat{p}}_{\text{positive}})}{n_{\text{trials}}}}).$$

For each ROC curve, we apply these two boundaries to calculate the shaded area as CI. ROC curves for Venus and Earth with different numbers of PC images to combine into one co-added image $N_{\mathrm{im}}$ (with different integration times) are shown in Fig. 7. As Venus is brighter than Earth, the performance for Venus is better than that for Earth. Increasing the integration time has the similar effect as increasing the planet intensity because both increase the expected number of photons arriving on the pixel. Moreover, the performance for both Venus and Earth is better with the higher number of PC images in the co-added image.

Fig. 7

ROC with CI for Venus and Earth using GLRT with different integration times. (a) ROC for co-added images each from $N_{\mathrm{im}}=200$ PC images. The curves for an integration time of 10 s overlaps with the plot’s boundary, which indicates perfect detection performance. ROC for co-added images each from (b) $N_{\mathrm{im}}=700$ PC images and (c) $N_{\mathrm{im}}=2000$ PC images. The shaded region behind each ROC curve is its 95% CI.

3.6.

Optimal Number of PC Images for One Co-Added Image

As we mentioned before, the number of PC images to combine into one co-added image $N_{\mathrm{im}}$ is an important hyperparameter to be chosen before doing detection via most of the methods. For example, GLRT’s performances vary with the number of PC image for one co-added image, as shown in Fig. 7. With too small $N_{\mathrm{im}}$, the TP rate and FA rate for the detection may not be desirable. With too big $N_{\mathrm{im}}$, precious observation time would be wasted. To choose the best $N_{\mathrm{im}}$, we can utilize the ROC curves.

First, given integration time, three parameters need to be specified: the minimum planet intensity to be detected, the maximum FA rate that can be accepted, and the minimum TP rate that is acceptable. Then for a different $N_{\mathrm{im}}$, the ROC is calculated via Monte Carlo simulation. Finally, the minimum $N_{\mathrm{im}}$ that can reach the requirements are chosen. For example, if we assume that the dimmest planet has the same intensity as Earth, the maximum acceptable FA rate is 0.16 and the minimum acceptable TP rate is 0.85. The acceptable FA rate and TP rate pairs are in the shaded green region in Fig. 8. We calculate the ROCs with different $N_{\mathrm{im}}$ for integration time 1 s and find that the ROC of $N_{\mathrm{im}}=700$ is the first one to reach the green region in Fig. 8. Thus we choose $N_{\mathrm{im}}=700$ as the optimal number of PC images to be co-added into the final image.

Fig. 8

ROC with CI for Earth using GLRT with different $N_{\mathrm{im}}$. The horizontal green dotted line is the acceptable TP rate and the vertical one is the acceptable FA rate region. Thus the acceptable true TP and acceptable FA region is the shaded green area. The ROC with $N_{\mathrm{im}}=700$ is the first ROC reaching the acceptable region.

3.7.

Comparison with Other Methods

We also compared GLRT with the performance of the detection method based on SNR map implemented in pyKLIP.³⁶ The algorithm computes the standard deviation in concentric annuli after masking the signal area in question as the level of noise in SNR. The width of the annuli used is the diameter of the PSF core in this paper. The SNR maps for the three example co-added images in Figs. 5(b) and 6(b) are shown in Fig. 9 to help visually compare the performance with the GLRT method. ROC curves are also calculated, shown in Fig. 10. The calculation uses the same set of images as the ones used for Fig. 7. As it is hard to visually compare the curves, we list the area under the ROC curve (AUC) for all the curves in Table. 3. AUC is an aggregate measure of performance across all possible thresholds. It can be interpreted as the probability that the model ranks a random positive example higher than a random negative example. AUC is 1 if the model’s decisions are all correct and 0 if all wrong. The disadvantage of this SNR definition is that the standard deviation calculation will be biased by the presence of point sources in annuli. In our case, the radius of Venus and Earth from the image center is close, so the signals are partially in each other’s annuli for the calculation of the noise standard deviation. Thus the SNR is biased, which is validated by the deteriorating performance for Earth when signals become stronger due to increased integration time or $N_{\mathrm{im}}$, shown in Figs. 10(b) and 10(c). Overall, GLRT outperforms SNR method on images that have not been post-processed; it is beyond the scope of this work to examine the effects of post-processing on each detection method.

Fig. 9

The SNR map calculated with the pyKLIP package:³⁷ result for the image in (a) Fig. 5(b) and (b) Fig. 6(b).

Fig. 10

ROC with CI for Venus and Earth using pyKLIP with different integration times. (a)–(c) The result where Earth and Venus are at the real locations as shown in Fig. 5(a). (d)–(f) The result where we eliminate the impact of planets for each other: when calculating the ROCs for Earth, we did not include Venus in the images; when calculating the ROCs for Venus, we did not include Earth in the images. (a), (d) ROC for co-added images each from $N_{\mathrm{im}}=200$ PC images; (b), (e) ROC from $N_{\mathrm{im}}=700$ PC images; and (c), (f) ROC from $N_{\mathrm{im}}=2000$ PC images.

Table 3

Comparison of AUC for GLRT and SNR method from pyKLIP.36 Unimpacted means only one planet is included in the image simulation.

	Venus 10 s	Earth 10 s	Venus 1 s	Earth 1 s	Venus 0.5 s	Earth 0.5 s
GLRT, 200 PC	1	1	0.9883	0.7374	0.8880	0.5797
GLRT, 700 PC	1	1	1	0.9503	0.9963	0.7490
GLRT, 2000 PC	1	1	1	0.9987	1	0.9275
SNR, 200 PC	1	0.8299	0.9714	0.6953	0.8804	0.6219
SNR, 700 PC	1	0.5059	0.9993	0.7159	0.9799	0.6863
SNR, 2000 PC	1	0.0325	1	0.4554	0.9985	0.6112
Unimpacted SNR, 200 PC	1	0.9999	0.9768	0.7595	0.8870	0.6529
Unimpacted SNR, 700 PC	1	1	0.9997	0.9069	0.9848	0.7711
Unimpacted SNR, 2000 PC	1	1	1	0.9806	0.9995	0.8753

4.

Iterative GLRT for Exozodiacal Dust

Exozodiacal dust is debris in the habitable zones of stars believed to come from extrasolar asteroids and comets.¹² Though the true structure of exozodiacal dust is unknown, as a first attempt, we simply assume it is axisymmetric, which is believed to be a reasonable approximation for small dust grains.³⁸ When the intensity of exozodiacal dust is similar to that of a target planet, the methods mentioned above have difficulty detecting the planet’s signal. We develop here an iterative GLRT to tackle this problem. It is essentially an expectation–maximization algorithm.³⁹ The planets’ signals and the exozodiacal dust are both unknown and need to be estimated. However, it is hard to estimate them accurately at the same time, and their estimation can influence each other. The solution is to iteratively estimate either the planets’ signal or the exozodiacal dust first and then use the estimation as a known factor to estimate the other until both estimates converge.

When there is exozodiacal dust, the model in Eq. (1), specified at pixel $(x,y)$, becomes

The exozodiacal dust degrades the detection and estimation. For example, if we directly apply the GLRT method introduced in the previous section for Fig. 4(d), we get a confusing FA rate map in Fig. 11(a); it is hard to distinguish Earth from the exozodiacal dust.

Fig. 11

Example of iterative GLRT applied to Fig. 4(d): (a) $T$ value map for Fig. 4(d). (b) The median value for each radius of Fig. 4(d), i.e., the exozodiacal dust estimation ${\mathit{d}}^{*}$ at initial step. (c) The residual ${\mathit{I}}_{\mathit{b}}$ after subtracting $d^{*}$ from $\mathit{I}$, i.e., the residual after subtracting (b) from Fig. 4(d). It is the initial estimation for the underlying image with only planets. (d) $T$ map for (c). (e) The new estimate of planets $\sum _{i=1}^{N_x}\sum _{j=1}^{N_y}[{\alpha }_{i,j}{\mathit{P}}_{i,j}(x,y)]$. After applying GLRT on (c) and get detection, we also obtained the intensity and position estimation of the planets. (f) Exozodiacal dust $I_p$ after subtracting estimated planet signals (e) and estimated local background from the original image Fig. 4(d). (g) The dust estimation at the final step. (h) The final residual ${\mathit{I}}_{\mathit{b}}$, i.e., the residual after subtracting (g) from Fig. 4(d). It is the final estimation for the underlying image with only planets and estimated local background. (i) $T$ value map for (h).

To reduce the number of unknowns, we assume that the dust is nearly axisymmetric, which may be a reasonable approximation for small dust grains.³⁸ Thus

Fig. 12

The flowchart describes the process of iterative GLRT. The blue box is the initialization step.

Table 4

Intensity and position estimation error for Fig. 4(f) via iterative GLRT methods

Fig. 13

The ROCs for the iterative GLRT. The ROCs without dust from Fig. 7 are added here for reference. The integration time here is 1 s. ROC for co-added images each from (a) $N_{\mathrm{im}}=200$ PC images; (b) $N_{\mathrm{im}}=700$ PC images; and (c) $N_{\mathrm{im}}=2000$ PC images.

5.

Conclusion

A starshade is a promising instrument for the direct imaging of Earth-like planets. In this paper, we briefly describe our process for simulating realistic starshade images and preliminary study of signal detection in starshade images, which no previous work has looked into. The core detection and estimation part are done by GLRT. We first obtain intensity estimates by QMLE. Then the likelihood ratio with respect to estimated parameters is calculated. After choosing an FA rate, we can threshold the image and detect the planets. For cases with exozodiacal dust, we split the process into two parts: dust estimation and signal estimation and use GLRT iteratively. Examples using these methods are shown in Secs. 3 and 4. The GLRT method successfully and efficiently flags potential planets with a concrete FA rate. It can help distinguish planet signals from artifacts caused by small starshade shape errors, such as a truncated petal tip. In addition, we provide a guidance to choose the best number of PC images to combine into one co-added image $N_{\mathrm{im}}$, utilizing the ROC curves. This will help utilize the observation time efficiently.

Due to the limitation of Gaussian approximation for the noise distribution in the image, Gaussian GLRT introduces detection performance improvement but not drastically, compared to the SNR method commonly used in high-contrast imaging. We have worked on an improved version of the GLRT method based on the accurate model for PC images rather than approximation and thus advances the detection performance;⁴⁰ we present the most recent result on this in Ref. 41. The performance can be further improved if we have prior knowledge about the probability distribution of the planets’ intensity in Eq. (21), which may be available after future exoplanet surveys. In this work, we present the iterative GLRT assuming face-on, uniform exozodiacal dust, but the same concept can be applied to more detailed models of the dust structure. In this work, spectral information is not discussed. However, the method introduced in this paper can be applicable to different cases. If images at different wavelengths are taken, the product of the likelihood at different wavelengths will be the final likelihood for these images. Then MLE and GLRT can be calculated, and thus a detection decision can be made.