Improving and Assessing Planet Sensitivity of the GPI Exoplanet Survey with a Forward Model Matched Filter

In this paper, we use 330 observations from the GPI Exoplanet Survey (GPIES) (Gemini programs GS-2014B-Q-500, GS-2015A-Q-500, and GS-2015B-Q-500; PI: B. Macintosh) to construct and test our matched filter. A typical GPIES epoch consists of 38  one-minute exposures in H band (1.5–1.8 μm). The number of usable raw spectral cubes can be lower because of degrading weather condition or isolated star tracking failures. We arbitrarily consider any data set with more than 20 usable exposures as complete. In some cases, a few exposures were added to the observing sequence in order to compensate for poor conditions, which resulted in a number of data sets with more than 38 exposures. For each star, we only consider the first complete epoch and ignore any follow-up observations that may have been made. We have not considered data sets with visible debris disks in order to avoid biasing the contrast curves. In H band, a GPI spectral cube has 37 wavelength channels and 281 × 281 pixels in the spatial dimensions, half of which are not filled with data, however, because of the tilted IFS field of view. Therefore, a typical GPIES data set includes approximately 1400 images at different position angles and wavelengths.

Spectral data cubes are built from raw IFS detector images using standard recipes from the GPI Data Reduction Pipeline³¹ version 1.3 and 1.4 (Perrin et al. 2016). The process includes correction for dark current, bad pixels, correlated read noise, and microphonics induced by cryocooler vibration (Ingraham et al. 2014). Flexure in the instrument slightly shifts the position of the lenslet micro-spectra on the detector. The offset is calibrated using argon arc lamp images at the target elevation, which are then compared with wavelength solution references taken at zenith (Wolff et al. 2014). Finally, each spectral cube is also corrected for optical distortion according to Konopacky et al. (2014). GPI images also contain four fainter copies of the unblocked PSF called satellite spots (Marois et al. 2006; Sivaramakrishnan & Oppenheimer 2006). The spots are used to estimate the location of the star behind the focal plane mask, and their flux allows the photometric calibration of the images. The GPI satellite spot-to-star flux ratio used here is $2.035\times {10}^{-4}$ for H band (Maire et al. 2014). In the following, the satellite spots are also median combined to estimate an empirical planet PSF that is wavelength dependent.

Simulated planets are necessary to optimize and characterize a detection algorithm because real point sources in high-contrast images are scarce. We decided to neglect for now the PSF smearing that is due to the sky rotation in a single exposure. The planetary spectra are selected from atmospheric models described in M. S. Marley et al. (2017, in preparation) and Saumon et al. (2012). For a given cloud coverage, the only model parameter with a significant effect on the shape of the spectrum in a single band is the temperature. However, this work is not about atmospheric characterization, so physically unrealistic model temperatures for a given object are not an issue as long as the shape of the spectrum matches. In the following, the references to T-type and L-type planets correspond to the analagous spectra of brown dwarfs. T-type spectra have strong methane absorption features, while L-type spectra are cloudy objects whose spectra are dominated by H₂O and CO. Of the known extrasolar planets, 51 Eridani b is an example of the T-type (Macintosh et al. 2015) and β Pictoris b of the L-type (Morzinski et al. 2015). The transmission spectrum of the Earth's atmosphere combined with that of the instrument is estimated by dividing the satellite spot spectrum with a stellar spectrum, which is then used to translate the spectrum from physical units into raw pixel values from the detector. The spectrum of the star is interpolated from the Pickles atlas, Pickles (1998), based on the star's spectral type.

Each individual image in the data set is speckle subtracted using a Python implementation of KLIP called PyKLIP³² (Wang et al. 2015). The KLIP algorithm consists of building and subtracting a model of the speckle pattern in an image, called science image, from a set of reference images that can be selected from the same data set or from completely different observations. In this paper, we use a combination of ADI and SDI strategies. The KLIP mathematical framework is summarized in Appendix A.1.1. A side effect of the speckle subtraction is the distortion of the planet PSF, referred to as self-subtraction. The self-subtraction is fully characterized in Appendix A.1.2.

All individual images are first high-pass filtered by subtracting a Gaussian-convolved image with a full-width at half-maximum (FWHM) of 12 pixels. Then, they are aligned and the reference images are scaled to the same wavelength as the science image. In order to account for spatial variation of the speckle behavior, KLIP is independently applied on small subsections of the field of view. Each image is therefore divided into small $100\ \mathrm{pixel}$ arcs to which a $10\ \mathrm{pixel}$ wide padding is added, as illustrated in Figure 1. For each sector, the reference library is built according to an exclusion criterion based on the displacement and the flux overlap (Marois et al. 2014) of the planet PSF between the science image and its reference images. The exclusion criterion and the reference library selection is described in Appendix A.1.3. The assumed spectrum for the companion used in the exclusion parameter is referred to as the reduction spectrum. In order to speed up the reduction, we only include the N_R = 150 most correlated images. This means that most of the images satisfying the exclusion criterion are in practice not used. In Section 4.1 it is shown that the exclusion criterion has a soft maximum around 0.7 for T-type planets, which is the value used in the following. The number of Karhunen-Loève modes of the reference library kept for the speckle subtraction is set to K = 30. This value has been chosen as a reasonable guess based on our experience, but it should be rigorously optimized in the future. In order to limit computation time, we have arbitrarily defined the outer working angle of the algorithm to $1^{\prime\prime}$ ($\approx 71$ pixels).

Zoom In Zoom Out Reset image size

In the field of signal processing, a matched filter is the linear filter maximizing the S/N of a known signal in the presence of additive noise (Kasdin & Braems 2006; Röver 2011; Cantalloube et al. 2015). A detailed description of the matched filter is presented in appendix Appendix A.2 and we only summarize the key results here. If the noise samples are independent and identically distributed, the matched filter corresponds to the cross-correlation of a template with the noisy data. In the context of high-contrast imaging, the pixels are neither independent nor identically distributed (i.e., heteroskedastic), which introduces a local noise normalization in the expression of the matched filter.

In a data set, each image is indexed by its exposure number τ and its wavelength λ. We define the vector ${{\boldsymbol{p}}}_{l}$, with $l=(\tau ,\lambda )$, as a specific speckle-subtracted image. Similarly, we define the matched-filter template ${{\boldsymbol{m}}}_{l}$ as the model of the planet signal in the corresponding processed image normalized such that it has the same broadband flux as the star. The whitening effect of the speckle subtraction allows one to assume uncorrelated residual noise, which significantly simplifies the matched filter. However, this assumption is not perfectly verified, and its consequence is discussed in Section 3.6. The maximum likelihood estimate of the planet contrast at separation ρ and position angle θ is then given by

$\begin{eqnarray}&&\tilde{\epsilon }(\rho ,\theta )=\frac{{\displaystyle \sum }_{l}{\,{\boldsymbol{p}}}_{l}^{\top }{{\boldsymbol{m}}}_{l}/{\sigma }_{l}^{2}}{{\displaystyle \sum }_{l}{\,{\boldsymbol{m}}}_{l}^{\top }{{\boldsymbol{m}}}_{l}/{\sigma }_{l}^{2},}\end{eqnarray} \tag{ 1 }$

where ${\sigma }_{l}$ is the local standard deviation at the position $(\rho ,\theta )$, assuming that it is constant in the neighborhood of the planet. Note that the planet model ${{\boldsymbol{m}}}_{l}$ approaches zero rapidly when moving away from its center $(\rho ,\theta )$, allowing one to only consider postage-stamp sized images containing the putative planet instead of the full images in Equation (1). Then, the theoretical S/N of the planet can be written

$\begin{eqnarray}&&{ \mathcal S }(\rho ,\theta )=\frac{{\displaystyle \sum }_{l}{\,{\boldsymbol{p}}}_{l}^{\top }{{\boldsymbol{m}}}_{l}/{\sigma }_{l}^{2}}{\sqrt{{\displaystyle \sum }_{l}{\,{\boldsymbol{m}}}_{l}^{\top }{{\boldsymbol{m}}}_{l}/{\sigma }_{l}^{2}}.}\end{eqnarray} \tag{ 2 }$

A detection can be claimed when the S/N is such that the observation cannot be explained by the null-hypothesis.

The calculation of the matched-filter template ${{\boldsymbol{m}}}_{l}$ is complicated by the distortion of the planet PSF after speckle subtraction. The characteristic effect of the self-subtraction manifests as negative wings on each side of the planet PSF and a narrowing of its central peak. It is common practice to define the matched-filter template as a 2D Gaussian, but this model fails to include most of the information about the distortion. In this paper, we present a novel matched-filter implementation using KLIP and the forward model of the planet PSF from Pueyo (2016) illustrated in Figure 2. The forward model is a linearized closed-form estimate of the distorted PSF, which is detailed in Appendix A.3. It is built from the same satellite spot-based PSF as the simulated planet injection in Section 2.3 and from the reduction spectrum. The forward model is a function of the PSF location on the image and therefore needs to be calculated for each pixel.

Zoom In Zoom Out Reset image size

In this paper, we define the throughput as the ratio of the estimated contrast of a planet after speckle subtraction over its original contrast. This definition might differ from the conventional idea of throughput in the data processing community for direct imaging. It is a measure of the ability of an algorithm to recover an unbiased contrast estimate of the point source. This can then be used to validate our implementation of the forward model. A low throughput indicates that the planet signal has been distorted during the speckle subtraction. The flux of each point source is estimated with Equation (1) and then compared to the known true contrast of the simulated planet. Figure 3 compares the throughput as a function the exclusion criterion, when using either the forward model or only the original PSF as a template. The estimation algorithm is otherwise identical in both cases. Figure 3 demonstrates that the forward model accurately models the self-subtraction because the throughput is close to one even for very aggressive reduction, while it drops very quickly to zero with the original PSF. However, as the exclusion parameter becomes too small, the linear approximation of the forward model breaks down and the throughput drops. As expected, the scatter in Figure 3 decreases with lower values of the exclusion criterion. The test was performed with two representative data sets at two different separations: (a) HR 7012 at $0\buildrel{\prime\prime}\over{.} 3$ and $0\buildrel{\prime\prime}\over{.} 6$, and (b) HD 131435 at $0\buildrel{\prime\prime}\over{.} 4$ and $0\buildrel{\prime\prime}\over{.} 8$. The data sets were chosen to represent different regimes of data quality as measured by their exoplanet contrast sensitivity in the fully reduced image; HR 7012 is an example of a good data set, while HD 131435 is of medium quality. A total of 35 simulated exoplanets were injected at each separation in seven copies of the same data set. For each copy of the data set, the five simulated planets, which are 72° apart, were rotated by 10° in position angle to better sample the image.

Zoom In Zoom Out Reset image size

The new algorithm is compared to two simpler matched filters, all using the same KLIP implementation, but involving a simple Gaussian PSF and no modeling of the planet self-subtraction. The three algorithms are referred to in this paper as Gaussian cross-correlation (GCC), Gaussian matched filter (GMF), and forward model matched filter (FMMF). The three methods are detailed in this section and illustrated in Figure 4. The different methods also differ in the way the data set is collapsed before the matched filter is performed. It is beyond the scope of this paper to evaluate the effect of each particular difference. All the algorithms presented in this paper are available in PyKLIP (Wang et al. 2015). A significant fraction of the code is shared with the Bayesian KLIP-FM Astrometry (BKA) method developed in Wang et al. (2016).

Zoom In Zoom Out Reset image size

The Gaussian cross-correlation (GCC) is the baseline algorithm because it is commonly used in high-contrast imaging. The overlapping sectors are mean combined after speckle subtraction in order to limit edge effects. Then, the processed data set is first derotated, coadded, and then collapsed using the reduction spectrum. The resulting image is cross-correlated with a 2D Gaussian kernel with a FWHM of 2.4 pixels. The size of the Gaussian has been chosen to be significantly smaller than the FWHM of the original PSF ($\approx 3.5\ \mathrm{pixels}$) to account for the self-subtraction. The cross-correlation is nothing more than a matched filter where the noise is spatially identically distributed. This assumption is not verified for GPI images, as discussed in Section 5.1. In practice, the noise properties only need to be azimuthally uniform; the theoretical S/N from the matched filter is always rescaled as a function of separation by estimating the standard deviation in concentric annuli, as discussed in Section 3.6. To summarize the Gaussian cross-correlation for GPIES, the template is a 20 × 20 pixel stamp projected on a 281 × 281 pixel image.

The Gaussian matched filter (GMF) projects the PSF template on the coadded cubes directly, as illustrated in Figure 4. Therefore, the template is a stamp cube with both spatial and wavelength dimensions. It also uses a Gaussian PSF with a 2.4 pixel FWHM, and the reduction spectrum is used to scale the template as a function of wavelength. The matched filter is calculated according to Equation (2) considering a single combined cube. As a consequence, this implementation does not assume identically distributed noise. The local noise standard deviation is estimated at each position in a 20 × 20 pixel stamp from which a central disk with a $2.5\ \mathrm{pixel}$ radius is removed. To summarize the GMF for GPIES, the template is a $20\times 20\times 37$ pixel cube stamp, which includes the spectral dimension, projected on a $281\times 281\times 37$ pixel datacube.

The forward model matched filter (FMMF) is performed on an uncollapsed data set according to Equation (2), as illustrated in Figure 4. It uses the forward model as a template. The sector padding provides the necessary margin to perform the projection of the template anywhere in the image. The local standard deviation in each image is estimated at the position $(\rho ,\theta )$ in a local arc defined by $[\rho -{\rm{\Delta }}\rho ,\rho +{\rm{\Delta }}\rho ]\times [\theta -{\rm{\Delta }}\rho /\rho ,\theta +{\rm{\Delta }}\rho /\rho ]$ and ${\rm{\Delta }}\rho =10\ \mathrm{px}$. In this case, the center of the arc is not masked, which results in an overestimation of the standard deviation in the presence of a planet signal. This does not significantly change the detectability of real objects, because the overestimated local standard deviation plays against both real objects and false positives.

In practice, only the values of the inner products ${{\boldsymbol{p}}}_{l}^{\top }{{\boldsymbol{m}}}_{l}$, ${{\boldsymbol{m}}}_{l}^{\top }{{\boldsymbol{m}}}_{l}$ and ${\sigma }_{l}(\rho ,\theta )$ are saved while the sectors are speckle subtracted in order to limit computer memory usage. The final sum of Equation (2) is performed at the very end. One advantage of not combining the data is that it is not necessary to derotate the speckle-subtracted images, as long as one accounts for the movement of the model in the data as a function of time and wavelength. This removes the image interpolation associated with derotation and therefore limits interpolation errors. The matched-filter calculation is run on a discretized grid, in this case, centered on each pixel in the final image. The matched filter is therefore slightly less sensitive to planets that are not also centered on a pixel. Assuming spatially randomly distributed planets, we estimate that the discretization results in an average loss of S/N of a few percentage points, compared to the case where every planet is centered on a pixel. To summarize the FMMF for GPIES, the template is a $20\times 20\times 37\times 38$ pixel multidimensional stamp, therefore it includes both spectral and time dimensions, projected on an uncombined $281\times 281\times 37\times 38$ pixel data set.

The FMMF reduction of a typical GPI 38 exposure data set requires around 30 wall clock hours on a computer equipped with 32 2.3 GHz cores and using ∼20 GB of random-access memory (RAM). As a consequence, a supercomputer is necessary to process an entire survey. Each data set is independent and can be run on separate nodes without sharing memory. In this paper, we have used the SLAC bullet cluster to process the entire campaign using on the order of 10⁶ CPU hours including the data processing necessary for the work presented in the remainder of the paper.

In practice, the theoretical S/N defined in Equation (2) needs to be empirically calibrated by estimating the standard deviation from the matched-filter map itself (Cantalloube et al. 2015). Indeed, the S/N is overestimated due to overly optimistic assumptions on the noise distribution. The residual noise is mostly white and Gaussian in areas that are not dominated by speckle noise, but both assumptions break down close to the mask because of speckle noise, causing the matched filter to lose some validity, although it remains relatively effective nonetheless. It is also hard to estimate the local noise accurately, because the noise properties vary from pixel to pixel, adding a layer of uncertainty on the theoretical S/N. This is why an empirical standard deviation is estimated in the matched-filter map as a function of separation using a 4-pixel wide annulus, as illustrated in Figure 5. In order to prevent a planet from biasing its own S/N, the standard deviation is calculated at each pixel while masking a disk with a 5-pixel radius centered on that pixel from the annulus. In addition, all the known astrophysical objects are also masked. In the particular case of injected simulated planets, the standard deviation is estimated from the planet free reduction. In the following, unless specified otherwise, any reference to an S/N relates to the calibrated S/N.

Zoom In Zoom Out Reset image size

For a centered Gaussian noise, the S/N is related to the false-positive rate (FPR) following

$\begin{eqnarray}&&\mathrm{FPR}=\displaystyle \frac{1}{2}\left(1-\mathrm{erf}\left(\displaystyle \frac{{\rm{S}}/{\rm{N}}}{\sqrt{2}}\right)\right).\end{eqnarray} \tag{ 3 }$

It is common practice to choose a $5\sigma$ detection threshold or higher, as discussed in Marois et al. (2008). For a Gaussian distribution, a $5\sigma$ threshold corresponds to an FPR equal to $\mathrm{FPR}=2.9\times {10}^{-7}$, which represents a false detection every $3.4\times {10}^{6}$ independent samples, or equivalently, on the order of 1000 GPI epochs. The deviation from Gaussianity of the residual noise in the final image will dramatically increase this false-positive fraction, as shown in Section 5.2. Small sample statistics (Mawet et al. 2014) also increases the false-positive fraction, but it is not expected to be a dominant term in this work and has therefore been neglected. Indeed, the inner working angle in the code is larger than three resolution elements at 1.6 μm.

Figure 6 gives an example of a reduction using the three different algorithms on the HR 7012 GPIES data set observed on 2015 April 08. This data set does not contain any visible astrophysical signal. Simulated planets were injected according to Section 2.3. With the exception of the anomaly at $0\buildrel{\prime\prime}\over{.} 4$, the FMMF consistently yields a higher S/N than the two other methods, and the final image shows fewer residual features. In the following, the simulated data sets are only reduced at the position of the simulated planets in order to speed up the algorithm. Figure 7 shows the processed data for several follow-up epochs of 51 Eridani also using the three algorithms. We show that FMMF would have marginally detected 51 Eridani b in all epochs, while this is not true for the GMF and the GCC.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The exclusion criterion is defined in Appendix A.1.3 and its optimal value is a trade-off between achieving a better speckle subtraction and maintaining a stronger planet signal. The forward model helps to keep the throughput close to unity for more aggressive reductions, therefore improving the overall S/N.

The optimal value of the exclusion criterion is found by calculating the S/N of simulated planets for different values of the parameter. We used similar simulated planets as for Figure 3 with a $1000\,{\rm{K}}$ cloud-free model spectrum and a $600\,{\rm{K}}$ cloud-free spectrum for the reduction. The choice of spectra is discussed in Section 4.2. Figure 8 shows that the exclusion parameter has a soft optimum around 0.7. The optimal exclusion criterion does not seem to significantly depend on separation or data set quality. This apparent stability of the FMMF optimum suggests that a single value of the exclusion criterion can be used for the entire survey. Consequently, all following reductions in H band will be performed with an exclusion criterion of 0.7. The optimal exclusion criterion is expected to change depending on the filter used for the observation, and future work should involve separate optimizations in the different spectral bands. The FMMF almost always yields a better S/N than the other two algorithms, but this does not necessarily mean that it has a better detection efficiency. Interestingly, neither GCC nor GMF have a consistent or even a well-defined optimum. It is very common in the field to reduce a data set with different sets of parameters and select the best one a posteriori, but we believe that the FMMF limits the need to fine-tune the parameters for each data set.

Zoom In Zoom Out Reset image size

In this section we discuss the optimization of the reduction spectrum used in the forward model and the reference library selection. The goal is to estimate the number of reduction spectra that should be used to recover the widest variety of planets. Currently detectable exoplanets and brown dwarfs are expected to feature spectra ranging from the T to the L spectral types. T-type objects are characterized by methane absorption bands visible in H band and an energy peak around 1.6 μm, while L-type objects feature a cloudy atmosphere with a flatter spectrum that peaks in the second half of the band.

We created ten copies of the HD 131435 GPIES data set, each with 16 simulated planets injected according to the spiral pattern from Figure 6. In each copy of the data set, the simulated planets were injected with one of ten spectra. The spectra were selected from a list of cloud-free and cloudy atmosphere models such that the most likely spectra from both T-type and L-type are represented. Note that the presence of clouds changes the temperature at which the methane features appear, but it does not significantly change the shape of the spectra. Then, each simulated data set was reduced with each of the ten reduction spectra, resulting in 100 different final products. For each of the planet spectra, the best reduction spectrum is defined by the one yielding the best median S/N for the 16 simulated planets. Figure 9 shows the median S/N of simulated planets as a function of their spectrum and the reduction spectrum. Surprisingly, the best reduction spectrum is not the one corresponding to the simulated spectrum. One possible explanation is that the same reduction spectrum is used for the forward model as for the reference library selection through the exclusion criterion. The spectrum is also a way to weight the spectral channels differently, which could effectively correct for a biased estimation of the standard deviation where there is planet signal. However, a deeper exploration of this effect is beyond the scope of this paper. We also conclude that only two spectra, the cloud-free $600\,{\rm{K}}$ and the cloudy $1300\,{\rm{K}}$, are necessary to allow the detection of most giant planets without a significant loss of S/N. This result is consistent with a similar study in Johnson-Groh (2017) using TLOCI (Marois et al. 2014). Although the $600\,{\rm{K}}$ spectrum is the optimal reduction spectrum for T-type objects, the methane-induced peak in H band is unrealistically sharp. It has indeed not yet been observed in a real directly imaged exoplanet. However, the reduction spectrum and the spectrum of the simulated planets can be different. In order to be more representative of the observations, we use a $1000\,{\rm{K}}$ cloud-free model spectrum similar to 51 Eridani b for the T-type injected planet, which has a softer methane-induced peak. The L-type simulated planets use the same $1300\,{\rm{K}}$ cloudy model spectrum as for the reduction.

Zoom In Zoom Out Reset image size

For the remainder of the paper, the T-type reduction refers to a cloud-free $600\,{\rm{K}}$ reduction spectrum, while the L-type reduction refers to a cloudy $1300\,{\rm{K}}$ reduction model.

By combining the 330 GPIES H -band observations after removing any known astrophysical point sources, we have been able to estimate the PDF of the residual noise up to an unprecedented precision in the fully reduced S/N maps. Figure 10 compares the ideal Gaussian PDF with the PDF of the different algorithms calculated from the normalized histograms of the pixel values of all the S/N maps. The statistics of the noise strongly deviates from the Gaussian distribution at occurrence rates much greater than the frequency of planets (Nielsen et al. 2013; Bowler et al. 2015; Galicher et al. 2016; Vigan et al. 2017), demonstrating the high occurrence of false positives with high S/N in direct imaging. The Gaussian cross-correlation has the same PDF as the S/N maps calculated from the speckle-subtracted images with no cross-correlation or matched filter. The GMF and the FMMF both significantly improve the statistics of the residuals but remain quite remote from an ideal Gaussian distribution. The excess of high S/N occurrences can be explained by either a truly non-Gaussian statistics or by a poor estimation of the standard deviation of the noise when calculating the S/N. For example, an underestimated standard deviation for a pixel will result in more high S/N false positives. The relative improvement between the cross-correlation and the matched filter suggests that the latter explanation may still be dominant, which we discuss in more detail in the next paragraph.

Zoom In Zoom Out Reset image size

Figure 11 shows the spatial densities of $3\sigma$ false positives as a function of separation and position angle. Owing to the correlation in the final S/N maps, one should count the number of speckles and not the raw number of pixels above a given threshold in order to count the number of false positives. Indeed, for any high S/N false detection, the size of the bump, and therefore the number of pixels above the threshold, will depend on the correlation length of the noise, which also depends on the algorithm used. The detection of false positives is therefore done recursively as follows. The highest S/N pixel is flagged, and a 4-pixel area radius is masked around it. This process is repeated until the only false positives left are below a predefined S/N threshold. Both the GCC and the GMF exhibit strong radial variations of their false-positive densities at the position of the sector boundaries, as seen in the top row plots of Figure 11. In addition, the GCC has a significant excess of false positives around 90° and 270° in position angle. This feature can be explained by the excess of speckle noise on both sides of the focal plane mask in the direction of the wind, which we refer to as the wind-butterfly. The pattern is still visible after combining the entire survey because the wind direction overwhelmingly favors the northeast on Cerro Pachon at the Gemini South telescope. The wind-butterfly breaks the assumption of azimuthally identically distributed noise, which we use when estimating the standard deviation in concentric annuli. The wind-butterfly explains why the probability density function of the GCC in Figure 10 is significantly higher than the other matched filters. The GMF does not suffer from this effect because the matched filter includes a normalization with respect to the local standard deviation estimated around each pixel. The FMMF features a similar PDF, meaning that a similar S/N detection threshold should yield the same number of false positives. The real performance of each algorithm is studied in the next sections. While a cross-correlation is a common planet-detection approach, our analysis suggests that it can be ill-suited if the noise varies azimuthally. One should instead use the expression for the matched filter from Equation (2). An alternative approach would be to vary the S/N threshold for each data set and as a function of position, but the lack of local independent samples to estimate a position dependent PDF at low false-positive rates makes this endeavor very challenging. Indeed, one needs to have probed the PDF of the noise at high S/N in order to evaluate the false-positive probability of any detection. For example, $5\sigma$ events are sufficiently rare that their occurrence rate can only be estimated from the data of a entire survey and not from individual images.

Zoom In Zoom Out Reset image size

The FMMF L-type reduction has significantly more false positives near the mask than a T-type reduction. This is likely because there is a higher density of speckles near the mask and the L-type spectrum is a better match to them than a sharper spectrum. We have also seen in Figure 10 that for the three algorithms, the L-type PDF has wider tails than the T-type PDF, suggesting that the number of L-type false positives will be higher at constant S/N.

In general, an improvement in the S/N does not guarantee a better detection efficiency because the false-positive rate could increase in the mean time. It is therefore important to compare the number of detected planets to the number of false positives. For example, Figure 8 showed that the Gaussian cross-correlation tends to have slightly higher S/N than the GMF but we will show in this section that the cross-correlation leads to more false positives. For this reason, receiver operating characteristic (ROC) have become increasingly popular in direct imaging to compare different algorithms (Caucci et al. 2007; Choquet et al. 2015; Pairet et al. 2016; Pueyo 2016; R. Jensen-Clem et al. 2017, in preparation). A ROC curve compares the false-positive fraction to the true positive fraction, i.e., completeness, as a function of the detection threshold. Alternatively, we have decided to replace the false-positive fraction by the number of false detections per epoch integrated over the entire image since it is easily translatable to survey efficiency and astrophysical background occurrence rate. We have also chosen to fix the planet contrast and assume a uniform planet distribution in separation and position angle. Although this approach does not address the dependence of the ROC curve on separation and planet contrast, it is sufficient to evaluate the relative performance of each algorithm. Multidimensional ROC curves could be used, but the contrast curves calculated in Section 6 already include most of the relevant information. The ROC curves for the three algorithms are shown in Figure 12. Each ROC curve has been built by injecting 16 simulated exoplanets, using either a $1000\,{\rm{K}}$ cloud-free T-type spectrum, reduced with a $600\,{\rm{K}}$ cloud-free model spectrum, or a $1300\,{\rm{K}}$ cloudy L-type spectrum, reduced with a the same cloudy spectrum. All planets were injected at a 4 × 10⁻⁶ contrast in 330 GPIES H -band data sets using the spiral pattern illustrated in Figure 6. The advantage of using a large number of data sets is to marginalize over the conditions in any one particular epoch. Figure 12 shows that FMMF yields a better completeness at any false-positive rate. In addition, the S/N threshold corresponding to a given false-positive rate is always higher in the case of the cross-correlation than for the matched filters because of the larger tail in the PDF. For example, it has a T-type false-positive rate that is roughly six times higher at $5\sigma$ than the two other algorithms.

Zoom In Zoom Out Reset image size

The detection threshold should be defined by the number of false positives that can reasonably be followed up during the survey. We set this false-positive rate at 0.05 per epoch corresponding to 30 false positives for the entire survey. The S/N threshold therefore depends on the algorithm, as shown in Figure 12.

A.1.1. Formalism

Karhunen-Loève image processing (KLIP) (Soummer et al. 2012) is one of the most widely used speckle subtraction frameworks. The general concept of speckle subtraction algorithms is to use a library of reference images from which a model of the speckle pattern is computed and then subtracted for each science image, effectively whitening the noise. KLIP calculates the principal components of this reference library and filters out the higher order modes that contain more noise and more planet signal. The reference library is a set of images spatially scaled to the same wavelength and containing random realization of the correlated speckles. We use a general observing strategy by combining ADI and SDI observations, sometimes referred to as ASDI, which is possible with an IFS like GPI. This choice of strategy implies that the astrophysical signal will be found both in the reference images and the science image. As a consequence, the algorithm is subject to self-subtraction, which is discussed in Appendix A.1.2. A detailed description of the mathematical formalism underlying the KLIP algorithm can be found in Savransky (2015), which we briefly summarize here.

The ${N}_{R}\times {N}_{\mathrm{pix}}$ matrix ${\boldsymbol{R}}$ is defined by concatenating the vectorized mean subtracted reference images r_j, which are vectors with ${N}_{\mathrm{pix}}$ elements, such that

$\begin{eqnarray}&&{\boldsymbol{R}}={[{{\boldsymbol{r}}}_{1},{{\boldsymbol{r}}}_{2},...,{{\boldsymbol{r}}}_{{N}_{R}}]}^{\top }.\end{eqnarray} \tag{ 6 }$

Then, the image-to-image sample covariance matrix describing the correlation of the images in the reference library is described, up to a proportional factor, by

$\begin{eqnarray}&&{\boldsymbol{C}}={\boldsymbol{R}}{{\boldsymbol{R}}}^{\top },\end{eqnarray} \tag{ 7 }$

with its eigenvectors and eigenvalues ${{\boldsymbol{v}}}_{1},{{\boldsymbol{v}}}_{2},...,{{\boldsymbol{v}}}_{{N}_{R}}$ and ${\mu }_{1},{\mu }_{2},...,{\mu }_{{N}_{R}}$, respectively. The Karhunen-Loève images form the optimal orthonormal basis that best represents any realization of the noise in a least-squares sense. They are defined as

$\begin{eqnarray}&&{{\boldsymbol{z}}}_{k}=\displaystyle \frac{1}{\sqrt{{\mu }_{k}}}\displaystyle \sum _{j=1}^{{N}_{R}}{{\boldsymbol{v}}}_{k}[j]{{\boldsymbol{r}}}_{j}.\end{eqnarray} \tag{ 8 }$

The images can be written in a more compact form by defining the matrices ${{\boldsymbol{V}}}_{K}={[{{\boldsymbol{v}}}_{1},{{\boldsymbol{v}}}_{2},...,{{\boldsymbol{v}}}_{K}]}^{\top }$, ${{\boldsymbol{Z}}}_{K}={[{{\boldsymbol{z}}}_{1},{{\boldsymbol{z}}}_{2},...,{{\boldsymbol{z}}}_{K}]}^{\top }$ and the diagonal matrix ${\boldsymbol{\Lambda }}=\mathrm{diag}{({\mu }_{1},{\mu }_{2},...,{\mu }_{K})}^{\top }$, with $K\leqslant {N}_{R}$ being the number of selected Karhunen-Loève images kept for the subtraction. Then, the collection of horizontally stacked Karhunen-Loève images ${{\boldsymbol{Z}}}_{K}={[{{\boldsymbol{z}}}_{1},{{\boldsymbol{z}}}_{2},...,{{\boldsymbol{z}}}_{K}]}^{\top }$ is equivalently written as

$\begin{eqnarray}&&{{\boldsymbol{Z}}}_{K}=\sqrt{{{\boldsymbol{\Lambda }}}^{-1}}{{\boldsymbol{V}}}_{K}{\boldsymbol{R}}.\end{eqnarray} \tag{ 9 }$

The speckle subtraction consists in subtracting the projection of the science image i onto the Karhunen-Loève basis from itself, which is given by

$\begin{eqnarray}&&{\boldsymbol{p}}={\boldsymbol{i}}-\displaystyle \sum _{k=1}^{K}\langle {\boldsymbol{i}}| {z}_{k}\rangle {z}_{k}={\boldsymbol{i}}-{{\boldsymbol{Z}}}_{K}^{\top }{{\boldsymbol{Z}}}_{K}{\boldsymbol{i}},\end{eqnarray} \tag{ 10 }$

where p refers to the processed image. The science and the reference images are often part of the same data set, and the science image refers to the image of a specific exposure and wavelength from which the speckle pattern is being subtracted. Typically, speckle subtraction algorithms are performed independently on small sectors of the image to account for spatial variations of the speckle noise, instead of using the full image, but the formalism is identical.

A.1.2. Self-subtraction

In the this section, we review the effect of an astrophysical signal as a small perturbation to Equation (10), following Pueyo (2016). In this section only, variables defined in Appendix A.1.1 are assumed to be planet free, and the hat indicates the perturbed variable, such that

$\begin{eqnarray}&&\forall j\in [1,{N}_{R}],\quad {\hat{{\boldsymbol{r}}}}_{j}={{\boldsymbol{r}}}_{j}+{\rm{\Delta }}{{\boldsymbol{r}}}_{j},\end{eqnarray} \tag{ 11 }$

with ${\rm{\Delta }}{{\boldsymbol{r}}}_{j}=\epsilon \,{{\boldsymbol{a}}}_{j}$, where ${{\boldsymbol{a}}}_{j}$ is the normalized planet signal. The planet signal is defined by its amplitude , a spectrum and the instrumental PSF before speckle subtraction. It is more convenient to normalize the planet signal to the same brightness as the host star such that is in practice the contrast of the planet. Similarly, we write ${\hat{{\boldsymbol{z}}}}_{j}={{\boldsymbol{z}}}_{j}+{\rm{\Delta }}{{\boldsymbol{z}}}_{j}$, ${\hat{{\boldsymbol{Z}}}}_{K}={{\boldsymbol{Z}}}_{K}+{\rm{\Delta }}{{\boldsymbol{Z}}}_{K}$, $\hat{{\boldsymbol{i}}}={\boldsymbol{i}}+{\rm{\Delta }}{\boldsymbol{i}}={\boldsymbol{i}}+\epsilon \,{\boldsymbol{a}}$, and $\hat{{\boldsymbol{p}}}={\boldsymbol{p}}+{\rm{\Delta }}{\boldsymbol{p}}\approx {\boldsymbol{p}}+\epsilon \,{\boldsymbol{m}}$, where ${\boldsymbol{m}}$ is the planet model that will be used to build the matched-filter template ${\boldsymbol{t}}$ introduced in Equation (2).

Applying KLIP to the perturbed images gives

$\begin{eqnarray}&&\hat{{\boldsymbol{p}}}=\hat{{\boldsymbol{i}}}-{\hat{{\boldsymbol{Z}}}}_{K}^{\top }{\hat{{\boldsymbol{Z}}}}_{K}\hat{{\boldsymbol{i}}},\\ &&{\boldsymbol{p}}+{\rm{\Delta }}{\boldsymbol{p}}=({\boldsymbol{i}}+\epsilon \,{\boldsymbol{a}})\\ &&\,-\,{({{\boldsymbol{Z}}}_{K}+{\rm{\Delta }}{{\boldsymbol{Z}}}_{K})}^{\top }\times ({{\boldsymbol{Z}}}_{K}+{\rm{\Delta }}{{\boldsymbol{Z}}}_{K})({\boldsymbol{i}}+\epsilon \,{\boldsymbol{a}}),\end{eqnarray} \tag{ 12 }$

and expanding Equation (12), one identifies ${\boldsymbol{p}}$, which should only contain residual speckles such that

$\begin{eqnarray}&&{\boldsymbol{p}}={\boldsymbol{i}}-{{\boldsymbol{Z}}}_{K}^{\top }{{\boldsymbol{Z}}}_{K}{\boldsymbol{i}}\approx 0.\end{eqnarray} \tag{ 13 }$

Then, the perturbed processed image is given by

$\begin{eqnarray}{\rm{\Delta }}{\boldsymbol{p}} & = & \epsilon \,{\boldsymbol{a}}\quad \mathrm{Original}\ \mathrm{planet}\ \mathrm{signal},\\ & & -\epsilon {{\boldsymbol{Z}}}_{K}^{\top }{{\boldsymbol{Z}}}_{K}{\boldsymbol{a}}\quad \mathrm{Oversubtraction},\\ & & -({{\boldsymbol{Z}}}_{K}^{\top }{\rm{\Delta }}{{\boldsymbol{Z}}}_{K}+{({{\boldsymbol{Z}}}_{K}^{\top }{\rm{\Delta }}{{\boldsymbol{Z}}}_{K})}^{\top }){\boldsymbol{i}}\quad \mathrm{Self} \mbox{-} \mathrm{subtraction},\\ & & +{ \mathcal O }({\epsilon }^{2})\quad \mathrm{Neglected}\ \mathrm{second} \mbox{-} \mathrm{order}\ \mathrm{terms}.\end{eqnarray} \tag{ 14 }$

The second and third terms of Equation (14) are the first-order distortions introduced by KLIP on the planet signal, as shown in Figure 2(b). The oversubtraction comes from the projection of the planet onto the unperturbed Karhunen-Loève modes, while the self-subtraction is the result of the projection of the speckles onto the perturbation of the modes. The former is unavoidable in a least-squares approach. The characteristic patterns for the self-subtraction are negative lobes in the radial and azimuthal direction that are due to the movement of the planet in the reference frames with wavelength and parallactic angle. The self-subtraction vanishes when there is no planet signal in the reference library as ${\rm{\Delta }}{{\boldsymbol{Z}}}_{K}=0$, which is the case when using a reference star differential imaging (RDI) approach. When the locally optimized combination of images (LOCI) of  Lafrenière et al. (2007a) is used, another way to avoid self-subtraction is to use optimization and subtraction region that do not overlap (see, Marois et al. 2010). However, the price of these strategies is a possibly less effective speckle subtraction.

A.1.3. Reference Library Exclusion Criterion

Part of the planet signal is lost in the processed image, especially when the planet PSF spatially overlaps in the reference images and in the science image. The distortion of the central peak of the PSF is more important when the overlap is large. A common practice to limit the distortion with ADI and SDI is to only include images in the reference library in which the planet has moved significantly compared to the science image. The importance of the distortion also depends on the relative brightness of the planet in the science image and in the references, which itself depends on the planet spectrum. Therefore, the selection of the reference images can be based on the relative flux overlap with the planet signal in the science image (Marois et al. 2014) instead of the sole consideration of the displacement.

The exclusion parameter chosen in this paper is a hybrid between a flux overlap and a displacement parameter. It corresponds to a pure displacement when comparing images at similar wavelength, but it more or less accepts images in the reference library when the wavelengths are different. This implementation of the exclusion criterion has given a significantly better result than any of the two other methods. The conversion between the three metrics is given for a specific planet PSF flux ratio between the science and the reference images in Figure 20. The selection of the reference images is illustrated in Figure 21 in a typical GPI data set. Only the 150 most correlated images with any given science frame are kept for the speckle subtraction in order to limit the size of the covariance matrix ${\boldsymbol{C}}$ and speed up the reduction. There is a stronger correlation between images with respect to wavelength than with respect to time, suggesting that SDI will give better results than ADI.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The choice of the exclusion criterion is a trade-off between an effective speckle subtraction and a limited self-subtraction. A low value of the exclusion criterion includes more correlated images in the reference library and results in a more aggressive reduction because of the higher degree of self-subtraction. The aggressiveness is also affected by the number of Karhunen-Loève modes used in the subtraction, but we keep this number fixed to K = 30 throughout this paper. We describe how the optimization of the exclusion parameter is performed in Section 4.1.

In the field of signal processing, a matched filter is the linear filter maximizing the S/N of a known signal in the presence of additive noise (Kasdin & Braems 2006). If the noise samples are independent and identically distributed, the matched filter corresponds to the cross-correlation of a template with the noisy data. In the context of high-contrast imaging, the pixels are neither independent nor identically distributed (i.e., heteroskedastic). Note that the matched filter is a very general tool that can be applied regardless of the algorithm used for speckle subtraction.

We define d as the vectorized speckle-subtracted data set representing a series of individual exposures of a star. Because GPI is an integral field spectrograph (IFS), the one-dimensional array d typically contains ${N}_{x}{N}_{y}{N}_{\lambda }{N}_{\tau }$ elements, where N_x and N_y are the number of pixels in the spatial dimensions (for a square image, N_x = N_y), N_λ is the number of spectral channels, and N_τ is the number of exposures. A typical GPIES data set is the concatenation of all the processed images $\{{p}_{l}\}$ from an epoch, with $l=(\tau ,\lambda )$, so ${N}_{x}{N}_{y}{N}_{\lambda }{N}_{\tau }=281\times 281\times 37\times 38$. If we assume that the data set contains a single point source, it can be decomposed as

$\begin{eqnarray}&&{\boldsymbol{d}}=\epsilon \,{\boldsymbol{t}}(\rho ,\theta )+{\boldsymbol{n}}.\end{eqnarray} \tag{ 15 }$

With ${\boldsymbol{t}}(\rho ,\theta )$ the vectorized template of the planet with separation and position angle $(\rho ,\theta )$. In practice we conveniently normalize the template to the same flux as the star such that is again the contrast. It is the concatenation of the planet model PSF for each image $\{{{\boldsymbol{m}}}_{l}\}$ and it has the same dimension as d. The second term n is for now assumed to be a Gaussian random vector assuming a zero mean and a covariance matrix ${\boldsymbol{\Sigma }}$. For conciseness, the $(\rho ,\theta )$ notation in ${\boldsymbol{t}}(\rho ,\theta )$ is dropped in the remainder of the paper, but the position dependence should be assumed for any point-source-related variables.

A matched filter can also be interpreted as a maximum likelihood ${ \mathcal L }(\epsilon ,\rho ,\theta )$ estimator of the planet amplitude and position assuming a Gaussian noise:

$\begin{eqnarray}&&{ \mathcal L }(\epsilon ,\rho ,\theta )=\displaystyle \frac{1}{\sqrt{2\pi | {\boldsymbol{\Sigma }}| }}\exp \left[\displaystyle \frac{-1}{2}{({\boldsymbol{d}}-\epsilon {\boldsymbol{t}})}^{\top }{{\boldsymbol{\Sigma }}}^{-1}({\boldsymbol{d}}-\epsilon \,{\boldsymbol{t}})\right].\end{eqnarray} \tag{ 16 }$

The most likely amplitude $\tilde{\epsilon }$ as a function of position is given by

$\begin{eqnarray}&&\tilde{\epsilon }(\rho ,\theta )=\frac{{{\boldsymbol{d}}}^{\top }{{\boldsymbol{\Sigma }}}^{-1}{\boldsymbol{t}}}{{{\boldsymbol{t}}}^{\top }{{\boldsymbol{\Sigma }}}^{-1}{\boldsymbol{t}}}.\end{eqnarray} \tag{ 17 }$

The theoretical S/N of a point source is given as a function of position by

$\begin{eqnarray}{ \mathcal S }(\rho ,\theta ) & = & \tilde{\epsilon }/{\sigma }_{\tilde{\epsilon }}\,= & \displaystyle \frac{{{\boldsymbol{d}}}^{\top }{{\boldsymbol{\Sigma }}}^{-1}{\boldsymbol{t}}}{\sqrt{{{\boldsymbol{t}}}^{\top }{{\boldsymbol{\Sigma }}}^{-1}{\boldsymbol{t}}}},\end{eqnarray} \tag{ 18 }$

which is a linear function of the square root of the maximized log-likelihood (Röver 2011; Cantalloube et al. 2015), where maximization is carried out over the point-source amplitude.

When we assume that the noise is uncorrelated and that its variance is constant in the neighborhood of the planet, the planet amplitude becomes

$\begin{eqnarray}&&\tilde{\epsilon }(\rho ,\theta )=\frac{{\displaystyle \sum }_{l}\,{{\boldsymbol{p}}}_{l}^{\top }{{\boldsymbol{m}}}_{l}/{\sigma }_{l}^{2}}{{\displaystyle \sum }_{l}{\,{\boldsymbol{m}}}_{l}^{\top }{{\boldsymbol{m}}}_{l}/{\sigma }_{l}^{2}},\end{eqnarray} \tag{ 19 }$

where ${\sigma }_{l}$ is the local standard deviation at the position $(\rho ,\theta )$. Note that the planet model ${{\boldsymbol{m}}}_{l}$ rapidly approaches rapidly when moving away from its center $(\rho ,\theta )$, allowing one to only consider postage-stamp-sized images containing the putative planet.

The theoretical S/N of the planet becomes

$\begin{eqnarray}&&{ \mathcal S }(\rho ,\theta )=\frac{{\displaystyle \sum }_{l}\,{{\boldsymbol{p}}}_{l}^{\top }{{\boldsymbol{m}}}_{l}/{\sigma }_{l}^{2}}{\sqrt{{\displaystyle \sum }_{l}{\,{\boldsymbol{m}}}_{l}^{\top }{{\boldsymbol{m}}}_{l}/{\sigma }_{l}^{2}}}.\end{eqnarray} \tag{ 20 }$

A detection can be claimed when the S/N is such that the observation cannot be explained by the null-hypothesis.

It is common practice to perform a matched filter with a simple Gaussian or a box template. Such templates are not optimal because they at best waste planet signal because of the self-subtraction. This paper uses the KLIP forward model from Pueyo (2016) to improve the matched-filter planet-detection sensitivity.

The perturbation of the processed image $\tfrac{{\rm{\Delta }}p}{\epsilon }$ from Equation (14) should be used as the template for the matched filter (also known as planet signature in Cantalloube et al. 2015 and in the forward model in Pueyo 2016). It can also be thought of as the derivative of the KLIP operator along the direction defined by the planet.

The ideal template $\tfrac{{\rm{\Delta }}p}{\epsilon }$ is a function of the original planet PSF and spectrum. The difficulty is the estimation of the Karhunen-Loève perturbation term ${\rm{\Delta }}{{\boldsymbol{Z}}}_{K}$ because it is a nonlinear function of the planet flux. Therefore Pueyo (2016) (E18) derived a first-order approximation of this term that can be used to calculate a close estimation of the planet PSF after speckle subtraction:

$\begin{eqnarray}\displaystyle \frac{{\rm{\Delta }}{{\boldsymbol{z}}}_{k}}{\epsilon } & = & -\left(\displaystyle \frac{{{\boldsymbol{v}}}_{k}^{\top }{{\boldsymbol{C}}}_{{AR}}{{\boldsymbol{v}}}_{k}}{2{\mu }_{k}}\right){{\boldsymbol{z}}}_{k}+\displaystyle \sum _{j=1,j\ne k}^{{N}_{R}}\left(\sqrt{\displaystyle \frac{{\mu }_{j}}{{\mu }_{k}}}\displaystyle \frac{{{\boldsymbol{v}}}_{j}^{\top }{{\boldsymbol{C}}}_{{AR}}{{\boldsymbol{v}}}_{k}}{{\mu }_{k}-{\mu }_{j}}\right){{\boldsymbol{z}}}_{j}\\ & & +\displaystyle \frac{1}{\sqrt{{\mu }_{k}}}{{\boldsymbol{A}}}^{\top }{{\boldsymbol{v}}}_{k}.\end{eqnarray} \tag{ 21 }$

With ${\boldsymbol{A}}={[{{\boldsymbol{a}}}_{1},{{\boldsymbol{a}}}_{2},...,{{\boldsymbol{a}}}_{{N}_{R}}]}^{\top }$ and ${{\boldsymbol{C}}}_{{AR}}={\boldsymbol{A}}{{\boldsymbol{R}}}^{\top }+{({\boldsymbol{A}}{{\boldsymbol{R}}}^{\top })}^{\top }$.

The forward model of the planet, also called template, is then given by

$\begin{eqnarray}&&{\boldsymbol{m}}={\boldsymbol{a}}-{{\boldsymbol{Z}}}_{K}^{\top }{{\boldsymbol{Z}}}_{K}{\boldsymbol{a}}-({{\boldsymbol{Z}}}_{K}^{\top }{\rm{\Delta }}{{\boldsymbol{Z}}}_{K}+{({{\boldsymbol{Z}}}_{K}^{\top }{\rm{\Delta }}{{\boldsymbol{Z}}}_{K})}^{\top })\displaystyle \frac{{\boldsymbol{i}}}{\epsilon },\end{eqnarray} \tag{ 22 }$

with ${\boldsymbol{m}}\approx \tfrac{{\rm{\Delta }}{\boldsymbol{p}}}{\epsilon }$. As a reminder, we have dropped the $(\rho ,\theta )$ upper script at the very beginning of this paper, but ${\boldsymbol{m}}$, ${\boldsymbol{a}}$ and ${\rm{\Delta }}{{\boldsymbol{Z}}}_{K}$ all depend on the position in the image.