Exoplanet Detection in Starshade Images

A starshade suppresses starlight by a factor of 1E11 in the image plane of a telescope, which is crucial for directly imaging Earth-like exoplanets. The state of the art in high contrast post-processing and signal detection methods were developed specifically for images taken with an internal coronagraph system and focus on the removal of quasi-static speckles. These methods are less useful for starshade images where such speckles are not present. This paper is dedicated to investigating signal processing methods tailored to work efficiently on starshade images. We describe a signal detection method, the generalized likelihood ratio test (GLRT), for starshade missions and look into three important problems. First, even with the light suppression provided by the starshade, rocky exoplanets are still difficult to detect in reflected light due to their absolute faintness. GLRT can successfully flag these dim planets. Moreover, GLRT provides estimates of the planets' positions and intensities and the theoretical false alarm rate of the detection. Second, small starshade shape errors, such as a truncated petal tip, can cause artifacts that are hard to distinguish from real planet signals; the detection method can help distinguish planet signals from such artifacts. The third direct imaging problem is that exozodiacal dust degrades detection performance. We develop an iterative generalized likelihood ratio test to mitigate the effect of dust on the image. In addition, we provide guidance on how to choose the number of photon counting images to combine into one co-added image before doing detection, which will help utilize the observation time efficiently. All the methods are demonstrated on realistic simulated images.


Introduction
A sun-like star is much brighter (typically 100 million to 10 billion times) than an Earth-like planet in its habitable zone. 1 Moreover, at a distance of 10 pc, the star and planets in its habitable zone are separated by around 0.1 arcseconds. 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, 2 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. 5 Thus, it imposes many challenging requirements on the telescope and instruments to design coronagraphs with ∼ 10 −10 starlight suppression; two mission concepts under study are Habitable Exoplanet Observatory (HabEx) 6 and Large UV/Optical/Infrared Surveyor (LUVOIR). 7 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). 8 In the mission, the starshade is launched separately and rendezvous with the telescope in orbit. Starshades are also baselined for the HabEx mission concept. 6 Starshades are a new technology, still in development. Coronagraphs, however, have been used on ground-based telescope 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 9 and unsuppressed starlight caused by errors in the starshade shape. 10 The scattered sunlight is confined to two extended lobes perpendicular to the direction of the sun and is constant during observations. 11 Its stability means 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. 12 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 post-processing (e.g., 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. Sec. 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.

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. 15  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), 17 which amplifies the signal in an electron-multiplication register. This process reduces the effective readout noise to less than one electron. 18 EMCCDs introduce an additional noise source: the multiplicative noise associated with the amplification process, which can be eliminated by 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 can not 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 less than one. In this way, photons aren't wasted. In this work, we choose the integration time for each PC image to be 1 second, 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 large enough dynamic range so the different intensity 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 im . In this work, we use N im = 2000, if not specified otherwise; guidance on how to choose N im is provided in Sec. 3.6.
The imaging field-of-view diameter is 9 arcseconds and each pixel is 0.021 arcsec by 0.021 arcsec.
To visualize the detector's performance, we use Monte Carlo to calculate the probability density functions (PDF) of the photon counts for different ground-truth photon flux 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 photon/s (our expected maximum rate) is well within the linear response regime.    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), 19 spectral diversity (SDI), 20 or multiple reference star images (RSDI) 21 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. Ref. 22  Recent studies also take the detection problem as a binary classification problem, using random-forest classifiers (SODIRF) and deep neural networks (SODINN). 30 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 generalized likelihood ratio test. 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.

GLRT model for the whole image
The model for an image containing multiple planet signals, background, and noise and is given by where I is the matrix for an image; α i,j is the intensity of signals at pixel (i, j) in unit of Jy (the value is zero if there is no signal at (i, j)); P i,j is the matrix denoting the values of a normalized PSF centered at (i, j); b is the matrix for background which contains star residual light and bias from the detector; ω is the matrix for noise; N x , N y are the number of pixels in x and y axis.
We assume that the noise in different pixels is an independent and identically distributed Gaussian random variable with mean zero and unknown constant variance. As each final image is the combination of many PC images, a Gaussian distribution should be a good approximation for the noise due to the central limit theorem. As a Gaussian distribution is used as an approximation for the true underlying distribution, the estimation is called Gaussian quasi-maximum likelihood estimation (QMLE). As long as the quasi-likelihood function is not oversimplified, the QMLE is consistent and asymptotically normal. It is less efficient than the MLE, but may only be slightly less efficient if the quasi-likelihood is constructed so as to minimize the loss of information relative to the actual likelihood. 31 We can then calculate the QMLE of α and b, which is equivalent to solving the optimization problem: This is an under-determined problem as we have 2N x N y unknown parameters, i.e., α, b and only N x N y data points, i.e., all the pixel values. Even assuming that we have a constant background in the whole image, we still have N x N y + 1 unknown parameters. Moreover, the assumption of constant background is not ideal as the background may have local features. To approach the problem, we use the idea of divide-and-conquer. If we assume no overlapping signals in the image, we are able to individually analyze smaller search areas that are the size of the PSF core.

GLRT model for a search area
In a small search area with the size of the PSF core, we can assume that there is only one planet signal. We also assume the background is constant, which should be a reasonable assumption over a small area. Moreover, noises are independent and identically distributed Gaussian random variables. Thus, the model for a search area is: where x i,j denotes the search area centered on the global pixel (i,j); b i,j is the background intensity; ω i,j denotes the Gaussian noise in this area. We only use the central core of the PSF for P i,j centered at (i,j). The model states that the area we are observing, x i,j , contains a signal centered at the center of this area along with constant background light and Gaussian noise. If there is no planet signal, α i,j should be zero. We stack pixel values in this target area into a column vector for easier mathematical manipulation. We do the same to the reference PSF and the noise matrix. The local model can be reformulated as a classical linear model: where x i,j is the vectorized target area centered at (i, j); N is the number of pixels in the area; b i,j is the constant background; P i,j (m) is the value at the m th pixel in the known vectorized template PSF centered at (i, j) for a point source signal of normalized intensity; α i,j is the planet's intensity, and w i,j is an N × 1 noise vector. θ i,j is unknown.
The conditional probability of this search area can be written as As the data x i,j are known and the parameters θ i,j , σ 2 i,j are unknown, this probability function is a likelihood function for the unknown parameters (so it is also denoted as L(θ i,j , σ 2 i,j ) above).
Taking the natural logarithm of both sides of Eq. (6), the log-likelihood of the search area in the co-added image is We maximize the log-likelihood function (equivalently maximizing the likelihood) to find the Maximum Likelihood estimator (the subscripts i, j for the estimated parameters are left out for simplicity): As mentioned previously, a Gaussian distribution is used as an approximation for the true underlying distribution; the estimation is also QMLE. The resulting estimation under H 1 is: And the estimated variance under H 1 is: Meanwhile, the QMLE under H 0 is: which is obtained from solving the constrained optimization problem: where A = [1,0]. And the estimated variance under H 0 is: In the problem of parameter estimation, we obtain information about the unknown parameters from the observed data of the random variables from the probability distribution governed by the parameters. The Fisher information matrix is a way to quantify the amount of information that the observable random variables carry about the unknown parameters. The definition of the Fisher information matrix is where the notation "Var" means the variance; "E" means expectation; φ is the vector of unknown parameters and φ = (α i,j b i,j σ 2 1 ) T for the alternative hypothesis H 1 . In our case (a linear model), the Fisher information matrix reduces to 32 Due to the block structure of I(φ), the variance-covariance ofθ G,1 can be estimated by I −1 and we can also derive the confidence intervals for the QMLE: 33 andb ± z (I −1 where z is the appropriate critical value (for example, 1.96 for 95 % confidence), and the notation ) ii means that we invert the matrix Iθ 1 first, and then take the ii component of the inverted matrix. The variance ofσ 2 1 is estimated by I −1 (σ 2 1 ) 32 where

Detection in a search area
The detection problem becomes a hypothesis testing problem: where A = [1,0].
To decide which hypothesis is true, it is intuitive to compare the posterior probability of H 1 and H 0 when x i,j occurs. We use Bayes' rule and define the odds ratio: where f (z) is the probability density function of the planet intensity and σ i,j is the standard deviation of the noise. We should decide H 1 is true if the odds ratio is larger than one. 22 However, we to a threshold. This test is called the generalized likelihood ratio test (GLRT). Moreover, not only is the planet intensity α i,j unknown, so are the background intensity and noise statistics. We use the QMLE mentioned above. Then, the ratio becomes whereθ i,j,k ,σ i,j,k are the estimated values under hypothesis H k , {k = 0, 1}.
To simplify further, H 1 is favored, i.e., a planet is more likely to exist at the test location, if where the threshold γ is based on the detection performance. 34 The probability of false alarms (false alarm rate) P F A and probability of detections P D are given by where Q is the probability of exceeding a given value; F 1,N −2 is an F distribution with one numerator degree of freedom and N-2 denominator degrees of freedom; and F 1,N −2 (λ i,j ) is a noncentral F distribution with one numerator degree of freedom and N-2 denominator degrees of freedom and noncentrality parameter λ i,j . 34 λ i,j is given by where θ i,j is the true value under H 1 and σ 2 is the true variance of the noise. This tells us that the probability of a false alarm only depends on the threshold, but the probability of detection depends on the planet intensity. The brighter the planet is, the higher the detection probability.

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 a false alarm 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.
When the pixel is at the boundary of the PSF in the image, 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 (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. The estimated planet intensity is the MLE of I at the estimated planet position. As we have shown in our previous work, 14 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 false alarm 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

Intensity Error = (Estimated Intensity−Real Intensity)
Real Intensity and report its value for all cases.

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 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 arcsec, so the position estimation is accurate to the pixel level. 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 im . Eq. (25) and Eq. (26) give the theoretical false alarm rate and true positive rate under a Gaussian assumption, and are therefore only an approximation. Thus, for a point (p f alse ,p positive ) on a ROC curve, we take the confidence interval using the two points: and  p f alse + z 1−α/2 For each ROC curve, we apply this two boundaries to calculate the shaded area as confidence interval. ROC curves for Venus and Earth with different numbers of PC images to combine into one co-added image N 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 more number of PC images in the co-added image.

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 im , is a important hyperparameter to be chosen before doing detection via most of the methods. For example, GLRT's performances varies with the number of PC image for one co-added image, as shown in Fig. 7. With too small N im , the true positive rate and false alarm rate for the detection may not be desirable. With too big N im , precious observation time would be wasted. To choose the best N 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 false alarm rate that can be accepted, and the minimum true positive rate that is acceptable. Then, for a different N im , the ROC is calculated via monte carlo simulation. Finally, the minimum N im that can reach the requirements are chosen. For example, if we assume the dimmest planet has the same intensity as Earth, the maximum acceptable false alarm rate is 0.16 and the minimum acceptable true positive rate is 0.85. The acceptable false alarm rate, true positive rate pairs are in the shaded green region in Fig. 8. We calculate the ROCs with different N im for integration time 1s and find that the ROC of N im = 700 is the first one to reach the green region in Fig. 8. Thus we choose N im = 700 as the optimal number of PC images to be co-added into the final image.

Comparison with other methods
We also compared GLRT with the performance of the detection method based on SNR map implemented in pyKLIP. 36 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 Fig. 5(b) and Fig. 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 im , shown in Fig. 10(b) and (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.

Iterative GLRT for Exozodiacal Dust
Exozodiacal dust is debris in the habitable zones of stars believed to come from extrasolar asteroids and comets. 12 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. 38 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  Fig. 5(b). (b) result for the image Fig. 6(b). GLRT to tackle this problem. It is essentially an Expectation-Maximization (EM) algorithm. 39 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: where d(x, y) is the exozodiacal dust at pixel (x, y) which is also unknown. We now have 3N x N y unknown parameters.
The exozodiacal dust degrades the detection and estimation. For example, if we directly apply the GLRT method introduced in previous section for Fig. 4(f), we get a confusing false alarm rate map in Fig. 11(a); it is hard to distinguish Earth from the exozodiacal dust.
To reduce the number of unknowns, we assume that the dust is nearly axisymmetric, which may be a reasonable approximation for small dust grains. 38 Thus, where r = x 2 + y 2 . Then Eq. (2) becomes We can not directly split the whole image into smaller areas and do detection separately like Eq. (9) as the estimation of dust in one area also depends on other areas at the same radii from the center.
To tackle this, we split the estimation of signals and exozodiacal dust into two steps. First, we take the median of the values for each radius to estimate the background. We use the median rather than mean to avoid the influence of the existence of planet signals at some radius, shown in Fig. 11(b).
This is equivalent to solving for the optimization problem for each r: The * denotes the estimate of the corresponding parameter. Then, we subtract the estimated background which contains bright exozodiacal dust, An example is shown in Fig. 11(c). Then, applying GLRT on this image I b produces the T value map in Fig. 11(d) and provides an estimation of the planets' positions and intensities, as shown in Fig. 11(e). The estimated planets are subtracted to get a better estimation of the background, as shown in Fig. 11(f) and the process is repeated iteratively. The procedure is summarized in Fig. 12 and the complete example is shown in Fig. 11. In Table 4, we summarize the intensity and position estimation error for the example. The ROC curves are shown in Fig. 13. The performance is undermined a little by the dust, compared to that without dust.

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 is done by GLRT. We first obtain intensity estimates by QMLE. Then, the likelihood ratio with respect to estimated parameters is calculated. After choosing a false alarm 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 Sec. 3 and Sec. 4. The GLRT method successfully and efficiently flags potential planets with a concrete false alarm 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 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; 40 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 detection decision can be made.

Disclosures
AH is a guest co-editor of this starshade special section.  images. The shaded region behind each ROC curve is its 95% confidence interval.
8 ROC with confidence interval for Earth using GLRT with different N im . The horizontal green dotted line is the acceptable true positive rate (TP) and the vertical one is the acceptable false alarm rate (FA) region. Thus, the acceptable true TP and acceptable FA region is the shaded green area. The ROC with N im = 700 is the first ROC reaching the acceptable region.

9
The SNR map calculated with the pyKLIP package 37 (a) result for the image in List of Tables   1 Parameters for simulation of solar system as viewed from 10 parsec.
2 Intensity estimation error and position estimation error comparison between results in images using different starshades.