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

TL;DR

This paper introduces a new matched filter algorithm with a forward model for direct exoplanet detection near bright stars, improving sensitivity and enabling comprehensive survey analysis.

Contribution

The paper develops a forward model matched filter incorporating PSF distortion, optimizing detection performance and enabling large-scale survey sensitivity assessment.

Findings

Detects 50% fainter objects than traditional methods at same false positive rate.

Uses only two spectral templates to recover any young Jovian exoplanets.

Assesses GPI-ES survey sensitivity as a function of star and companion spectral type.

Abstract

We present a new matched filter algorithm for direct detection of point sources in the immediate vicinity of bright stars. The stellar Point Spread Function (PSF) is first subtracted using a Karhunen-Lo\'eve Image Processing (KLIP) algorithm with Angular and Spectral Differential Imaging (ADI and SDI). The KLIP-induced distortion of the astrophysical signal is included in the matched filter template by computing a forward model of the PSF at every position in the image. To optimize the performance of the algorithm, we conduct extensive planet injection and recovery tests and tune the exoplanet spectra template and KLIP reduction aggressiveness to maximize the Signal-to-Noise Ratio (SNR) of the recovered planets. We show that only two spectral templates are necessary to recover any young Jovian exoplanets with minimal SNR loss. We also developed a complete pipeline for the automated detection of point source candidates, the calculation of Receiver Operating Characteristics (ROC), false positives based contrast curves, and completeness contours. We process in a uniform manner more than 330 datasets from the Gemini Planet Imager Exoplanet Survey (GPIES) and assess GPI typical sensitivity as a function of the star and the hypothetical companion spectral type. This work allows for the first time a comparison of different detection algorithms at a survey scale accounting for both planet completeness and false positive rate. We show that the new forward model matched filter allows the detection of $50\%$ fainter objects than a conventional cross-correlation technique with a Gaussian PSF template for the same false positive rate.

Tables

Table 1. Table 1: Definition of Scalar Variables.

Symbol	Type	Description
$N_x,N_y$	$\mathbb{N}\times \mathbb{N}$	Dimensions of the image.
$N_{\mathrm{pix}}$	$\mathbb{N}$	Number of pixels in a sector. If the sector is the entire image, then $N_{\mathrm{pix}}=N_x{N}_y$.
$N_{\lambda }$	$\mathbb{N}$	Number of wavelength channels in the IFS data cubes.
$N_{\tau }$	$\mathbb{N}$	Number of exposures in an epoch.
$N_R$	$\mathbb{N}$	Number of images in the reference library.
$K$	$\mathbb{N}$	Number of selected Karhunen-Loève images to be used in the speckle subtraction.
$\rho$	${ℝ}^{+}$	Separation to the star.
$\theta$	$ℝ$	Position Angle with the origin at the image North.
${\sigma }_l$	${ℝ}^{+}$	Local standard deviation of the noise at the position $(\rho ,\theta )$ in the speckle subtracted images.
$𝒮(\rho ,\theta )$	$ℝ$	Theoretical Signal-to-Noise Ratio (SNR) at the position $(\rho ,\theta )$ assuming independent Gaussian Noise. This is not the true SNR of the signal.
$ϵ$	${ℝ}^{+}$	Contrast of the point source relative to the host star. Different units could be used depending on the normalization of the planet model.
$\widetilde{ϵ}(\rho ,\theta )$ or $\widetilde{ϵ}$	$ℝ$	Maximum likelihood estimation of the amplitude of the planet signal $ϵ$ at the position $(\rho ,\theta )$.
${\sigma }_{ϵ}$	${ℝ}^{+}$	Standard deviation of the planet contrast at the position $(\rho ,\theta )$. ${\sigma }_{ϵ1}$ and ${\sigma }_{ϵ2}$ respectively refer to the first and follow-up observation of a star.
${\mu }_j$	$ℝ$	Eigenvalues of the matrix $𝑪$.
$\gamma$	$ℝ$	Conversion factor from the contrast matched filter values $𝒮$ to the contrast $\zeta$ such that $\zeta =\gamma 𝒮$.
$\eta$	$ℝ$	SNR detection threshold. ${\eta }_2$ is the detection threshold in the follow-up observation, which can be smaller.
${\sigma }_{\gamma }$	${ℝ}^{+}$	Standard deviation of $n_{\gamma }$
${\sigma }_{𝒮}$	${ℝ}^{+}$	Standard deviation of the noise in the matched filter map as a function of position.

Table 2. Table 2: Definition of Vectors.

Symbol	Length	Description
${𝐫}_j$	$N_{\mathrm{pix}}$	Vectorized images or sectors from the KLIP reference library.
${\widehat{𝐫}}_j$	$N_{\mathrm{pix}}$	Vectorized perturbed images or sectors from the KLIP reference library including some planet signal.
$\Delta {𝐫}_j$	$N_{\mathrm{pix}}$	Vectorized perturbation of the images or sectors from the KLIP reference library.
${𝐯}_j$	$N_R$	Eigenvectors of the matrix $𝑪$.
${𝐳}_k$	$N_{\mathrm{pix}}$	Karhunen-Loève images or sectors for the reference library matrix $𝑹$.
${\widehat{𝐳}}_k$	$N_{\mathrm{pix}}$	Perturbed Karhunen-Loève images or sectors for the reference library matrix $𝑹$ including planet signal.
$\Delta {𝐳}_k$	$N_{\mathrm{pix}}$	Perturbation of the Karhunen-Loève images or sectors for the perturbed reference library.
$𝐢$	$N_{\mathrm{pix}}$	Vectorized science image or sector before speckle subtraction.
$\widehat{𝐢}$	$N_{\mathrm{pix}}$	Vectorized perturbed science image or sector before speckle subtraction including planet signal.
$\Delta 𝐢$	$N_{\mathrm{pix}}$	Vectorized perturbation of the science image or sector before speckle subtraction.
${𝐚}_j$	$N_{\mathrm{pix}}$	Vectorized planet signal in the reference image or sector ${\hat{r}}_j$ before speckle subtraction.
$𝐚$	$N_{\mathrm{pix}}$	Vectorized planet signal in the science image or sector before speckle subtraction.
$𝐩$	$N_{\mathrm{pix}}$	Vectorized processed image or sector after speckle subtraction.
$\widehat{𝐩}$	$N_{\mathrm{pix}}$	Vectorized perturbed processed image or sector after speckle subtraction including planet signal.
$\Delta 𝐩$	$N_{\mathrm{pix}}$	Vectorized perturbation of the processed image or sector after speckle subtraction.
$𝐦$	$N_{\mathrm{pix}}$	Vectorized planet model in the processed image or sector after speckle subtraction.
$𝐝$	$N_x{N}_y{N}_{\lambda }{N}_{\tau }$	Vectorized speckle subtracted dataset.
$𝐭(\rho ,\theta )$ or $𝐭$	$N_x{N}_y{N}_{\lambda }{N}_{\tau }$	Vectorized planet template for the speckle subtracted dataset.
$𝐧$	$N_x{N}_y{N}_{\lambda }{N}_{\tau }$	Multivariate noise distribution of the speckle subtracted dataset. The noise is assumed to be Gaussian for any mathematical derivation.
${𝐩}_l$	$N_x{N}_y$	Vectorized speckle subtracted image with exposure $\tau$ and wavelength $\lambda$.
${𝐦}_l(\rho ,\theta )$ or ${𝐦}_l$	$N_x{N}_y$	Vectorized planet template for the speckle subtracted image with exposure $\tau$ and wavelength $\lambda$.

Table 3. Table 3: Definition of Matrices.

Symbol	Dimensions	Description
$𝚺$	${\left(N_x{N}_y{N}_{\lambda }{N}_{\tau }\right)}^2$	Covariance matrix of the noise $n$ in the speckle subtracted
$𝑹$	$N_R\times N_{\mathrm{pix}}$	Reference library matrix such that $𝑹={[r_1,r_2,\dots ,r_{N_R}]}^{\top }$.
$𝑪$	$N_R\times N_R$	Matrix proportional to the time covariance matrix defined as $𝑪=𝑹{𝑹}^{\top }$.
${𝑽}_K$	$K\times N_R$	Eigenvectors matrix defined as ${𝑽}_K={[v_1,v_2,\dots ,v_K]}^{\top }$.
${𝒁}_K$	$K\times N_{\mathrm{pix}}$	Karhunen-Loève images matrix defined as ${𝒁}_K={[z_1,z_2,\dots ,z_K]}^{\top }$.
${\widehat{𝒁}}_K$	$K\times N_{\mathrm{pix}}$	Perturbed Karhunen-Loève images matrix defined as ${\widehat{𝒁}}_K={[{\hat{z}}_1,{\hat{z}}_2,\dots ,{\hat{z}}_K]}^{\top }$ including planet signal.
$\Delta {𝒁}_K$	$K\times N_{\mathrm{pix}}$	Perturbation of the Karhunen-Loève images matrix defined as $\Delta {𝒁}_K={[\Delta z_1,\Delta z_2,\dots ,\Delta z_K]}^{\top }$.
$𝚲$	$K\times N_{\mathrm{pix}}$	Diagonal matrix defined as $𝚲=\text{diag}{({\lambda }_1,{\lambda }_2,\dots ,{\lambda }_{N_R})}^{\top }$.
$𝑨$	$N_R\times N_{\mathrm{pix}}$	Planet signal component in the reference library $𝑹$ and $𝑨={[a_1,a_2,\dots ,a_{N_R}]}^{\top }$.
${𝑪}_{𝑨𝑹}$	$N_R\times N_R$	Matrix used in the forward calculation defined as ${𝑪}_{𝑨𝑹}=𝑨{𝑹}^{\top }+{\left(𝑨{𝑹}^{\top }\right)}^{\top }$.

Equations

$$\tilde{\epsilon}(\rho,\theta)=\sum\limits_{l}{\bm{p}}_{l}^{\top}{\bm{m}}_{l}/\sigma_{l}^{2}\Bigg{/}\sum\limits_{l}{\bm{m}}_{l}^{\top}{\bm{m}}_{l}/\sigma_{l}^{2},$$

$$\mathcal{S}(\rho,\theta)=\sum\limits_{l}{\bm{p}}_{l}^{\top}{\bm{m}}_{l}/\sigma_{l}^{2}\Bigg{/}\sqrt{\sum\limits_{l}{\bm{m}}_{l}^{\top}{\bm{m}}_{l}/\sigma_{l}^{2}}.$$

$$\mathrm{FPR}=\frac{1}{2}\left(1-\mathrm{erf}\left(\frac{\mathrm{SNR}}{\sqrt{2}}\right)\right).$$

$$P_{\mathrm{detec}}(\mathrm{SNR},{\sigma }_{ϵ1}/{\sigma }_{ϵ2})={\int }_{\mathrm{s}}(1-\mathrm{CDF}({\eta }_2-s{\sigma }_{ϵ1}/{\sigma }_{ϵ2}))\mathrm{PDF}(s-\mathrm{SNR})\mathrm{d}s,$$

$${\sigma }_{ϵ}(ϵ)=\sqrt{{\gamma }^2{\sigma }_{\mathcal{S}}^2+\frac{{ϵ}^2}{{\gamma }^2}{\sigma }_{\gamma }^2},$$

$$\mathbf{R}={\left[{\mathbf{r}}_1,{\mathbf{r}}_2,\dots ,{\mathbf{r}}_{N_R}\right]}^{\top }.$$

$$\mathbf{C}=\mathbf{R}{\mathbf{R}}^{\top },$$

$${\mathbf{z}}_k=\frac{1}{\sqrt{{\mu }_k}}\sum _{j=1}^{N_R}{\mathbf{v}}_k[j]{\mathbf{r}}_j.$$

$${\mathbf{Z}}_K=\sqrt{{\mathbf{\Lambda }}^{-1}}{\mathbf{V}}_K\mathbf{R}.$$

$$\mathbf{p}=\mathbf{i}-\sum _{k=1}^K⟨\mathbf{i}|z_k⟩z_k=\mathbf{i}-{\mathbf{Z}}_K^{\top }{\mathbf{Z}}_K\mathbf{i},$$

$$\forall j\in [1,N_R],{\hat{\mathbf{r}}}_j={\mathbf{r}}_j+\Delta {\mathbf{r}}_j,$$

$$\hat{\mathbf{p}}$$

$$\mathbf{p}+\Delta \mathbf{p}$$

$$\mathbf{p}=\mathbf{i}-{\mathbf{Z}}_K^{\top }{\mathbf{Z}}_K\mathbf{i}\approx 0.$$

$$\Delta \mathbf{p}=$$

$$-$$

$$+$$

$$\mathbf{d}=ϵ\mathbf{t}(\rho ,\theta )+\mathbf{n},$$

$$\mathcal{L}(ϵ,\rho ,\theta )=\frac{1}{\sqrt{2\pi |\mathbf{\Sigma }|}}\exp \left[\frac{-1}{2}(\mathbf{d}-ϵ\mathbf{t}{)}^{\top }{\mathbf{\Sigma }}^{-1}(\mathbf{d}-ϵ\mathbf{t})\right].$$

$$\tilde{\epsilon}(\rho,\theta)={\bm{d}}^{\top}{\bm{\Sigma}}^{-1}{\bm{t}}\Big{/}{\bm{t}}^{\top}{\bm{\Sigma}}^{-1}{\bm{t}}.$$

$$\mathcal{S}(\rho ,\theta )$$

$$\displaystyle={\bm{d}}^{\top}{\bm{\Sigma}}^{-1}{\bm{t}}\Big{/}\sqrt{{\bm{t}}^{\top}{\bm{\Sigma}}^{-1}{\bm{t}}},$$

$$\tilde{\epsilon}(\rho,\theta)=\sum\limits_{l}{\bm{p}}_{l}^{\top}{\bm{m}}_{l}/\sigma_{l}^{2}\Bigg{/}\sum\limits_{l}{\bm{m}}_{l}^{\top}{\bm{m}}_{l}/\sigma_{l}^{2},$$

$$\mathcal{S}(\rho,\theta)=\sum\limits_{l}{\bm{p}}_{l}^{\top}{\bm{m}}_{l}/\sigma_{l}^{2}\Bigg{/}\sqrt{\sum\limits_{l}{\bm{m}}_{l}^{\top}{\bm{m}}_{l}/\sigma_{l}^{2}}.$$

$$\frac{\Delta {\mathbf{z}}_k}{ϵ}=-\left(\frac{{\mathbf{v}}_k^{\top }{\mathbf{C}}_{\mathbf{A}\mathbf{R}}{\mathbf{v}}_k}{2{\mu }_k}\right){\mathbf{z}}_k+\sum _{j=1,j\ne k}^{N_R}\left(\sqrt{\frac{{\mu }_j}{{\mu }_k}}\frac{{\mathbf{v}}_j^{\top }{\mathbf{C}}_{\mathbf{A}\mathbf{R}}{\mathbf{v}}_k}{{\mu }_k-{\mu }_j}\right){\mathbf{z}}_j+\frac{1}{\sqrt{{\mu }_k}}{\mathbf{A}}^{\top }{\mathbf{v}}_k,$$

$$\mathbf{m}=\mathbf{a}-{\mathbf{Z}}_K^{\top }{\mathbf{Z}}_K\mathbf{a}-\left({\mathbf{Z}}_K^{\top }\Delta {\mathbf{Z}}_K+{\left({\mathbf{Z}}_K^{\top }\Delta {\mathbf{Z}}_K\right)}^{\top }\right)\frac{\mathbf{i}}{ϵ},$$

Full text

instrumentation: adaptive optics — methods: statistical — planetary systems — surveys — techniques: high angular resolution — techniques: image processing.

††facilities: Gemini:South(GPI)††software: pyKLIP444Documentation available at http://pyklip.readthedocs.io/en/latest/. (Wang et al., 2015), GPI Data Reduction Pipeline555Documentation available at http://docs.planetimager.org/pipeline/. (Perrin et al., 2016), astropy666http://www.astropy.org (Astropy Collaboration et al., 2013)

1 Introduction

Direct imaging techniques spatially resolve exoplanets from their host star by using high-contrast imaging instruments usually combined with the power of large telescopes, adaptive optics, coronagraphs and sophisticated data processing. This technique currently allows the detection of young ($<300\,$Myr), massive ($>2\,M_{\mathrm{Jup}}$), self-luminous exoplanets at host-star separations not yet covered by indirect methods ($a>5\,$au) and therefore helps to constrain planet population statistics. The Gemini Planet Imager (GPI) (Macintosh et al., 2014), operating on the Gemini South telescope, is one of the latest generation of high contrast instruments with extreme adaptive optics. The GPI Exoplanet Survey (GPIES) is targeting 600 young stars and has to date observed more than half of them. As part of the survey, it imaged several known systems and discovered the exoplanet 51 Eridani b (Macintosh et al., 2015).

High contrast images suffer from spatially correlated noise, called speckles, which originate from optical aberrations in the instrument as well as a diffuse light component resulting from the time averaged uncorrected atmospheric turbulence. The correlation length of the speckles in a raw image is equal to the size of the unocculted Point Spread Function (PSF). Speckles are often described as quasi-static, since they are correlated across the observed spectral band and across the exposures of a typical observation (Perrin et al., 2003; Bloemhof et al., 2001; Sivaramakrishnan et al., 2002). It has been shown that the Probability Density Function (PDF) of the speckle noise is not Gaussian, but rather better described by a modified Rician distribution (Soummer & Aime, 2004; Bloemhof, 2004; Fitzgerald & Graham, 2006; Soummer et al., 2007; Hinkley et al., 2007; Marois et al., 2008; Mawet et al., 2014). The Rician PDF has a larger positive tail than a Gaussian distribution yielding a comparatively higher number of false positives at constant Signal-to-Noise Ratio (SNR).

In a single image, disentangling the signal of a faint planet buried under speckle noise is challenging because the spatial scale of speckle noise and the signal from a planet are similar, corresponding to the size of the Point Spread Function (PSF). However, speckles and astrophysical signals behave differently with time and wavelength. This diversity can be used to build a model of the speckle pattern and subtract it from all images. Two observing strategies are commonly used for speckle subtraction with instruments like GPI; Angular Differential Imaging (ADI) (Marois et al., 2006) and Spectral Differential Imaging (SDI) (Marois et al., 2000; Sparks & Ford, 2002). The observing setup for ADI is different from traditional imaging with altitude/azimuth telescopes, as the instrument field derotator is switched off or adjusted to keep the telescope pupil fixed with respect to sky rotation. As a consequence, the astrophysical signal rotates on the detector, following the sky rotation with respect to the telescope (i.e. parallactic angle evolution), while the speckle pattern remains relatively stable. In a similar tactic, SDI exploits the radial linear dependence of the speckle pattern with wavelength to separate it from the planet signal, which remains at the same position at all wavelengths. By definition, an observation with an Integral Field Spectrograph (IFS), like GPI, provides both temporal and spectral diversity of the speckle pattern for ADI or SDI to be used separately or in tandem.

The most common spectral subtraction algorithms are Locally Optimized Combination of Images (LOCI) Lafrenière et al. (2007a), which uses a least square approach to optimally subtract the speckle noise, and Karhunen-Loéve Image Processing (KLIP) Soummer et al. (2012), which regularizes the least square problem by filtering out the high order singular modes. Following speckle subtraction, point sources can be searched for by using a matched filter or a more general Bayesian model comparison framework as shown in Kasdin & Braems (2006). However, the distortion of the planet PSF caused by the speckle subtraction algorithm, referred to as self-subtraction, can make it difficult to define an accurate matched filter template. The self-subtraction can be accurately modeled in simplified cases, as shown in the matched filter approach by Cantalloube et al. (2015) using ADI subtracted pair images. In the context of planet characterization, the self-subtraction also biases the photometry and the astrometry of the object. The inverse problem is usually solved by injecting negative planets in the raw images and iteratively minimizing the image residuals (Marois et al., 2010; Morzinski et al., 2015). Recently, Pueyo (2016) derived a closed-form approximation of the self-subtraction in KLIP but without applying it in the context of a matched filter. Wang et al. (2016) used this new forward model in a Bayesian framework to estimate the astrometry of $\beta$ Pictoris b.

The main challenge with uniformly and systematically characterizing the detections of a high-contrast imaging exoplanet survey is the high number of false positives even at relatively high SNR. The detection threshold is hard to define because the PDF of the residual noise is generally unknown and depends on the instrument, the choice of data processing, and the dataset itself. The lack of well-known false positive rates makes it very difficult to evaluate the performance of algorithms relative to one another. Currently, candidates are discarded as false positives by visual inspection, which does not permit a rigorous calculation of planet completeness. In order to accurately characterize planet detection statistics, and ultimately constrain the underlying planet population in a uniform and unbiased manner, it is important to improve and systematize exoplanet detection methods.

The goal of this paper is to define a systematic and rigorous approach for exoplanet detection in large direct imaging surveys such that the long-period exoplanet population can be inferred in a meaningful manner. Using the latest KLIP framework, we develop an automated matched-filter based detection algorithm that includes a forward model of the planet self-subtraction (e.g., Pueyo, 2016) and accounts for the noise variations in the spatial, temporal and spectral dimensions of a dataset. As of the end of 2016, GPIES has already observed 330 stars, which allows us to precisely estimate the false positive rate and define meaningful detection thresholds for the entire survey. We conduct a rigorous set of tests to characterize state-of-the-art detection algorithms and demonstrate that a Forward Model Matched Filter (FMMF) most effectively recovers a planet signal while reducing the number of false positive detections. The paper is structured as follow:

• •

GPIES observations and data reduction are presented in Section 2,

• •

The matched filter is described in Section 3,

• •

The optimization of the reduction parameters for GPIES is presented in Section 4,

• •

The residual noise is characterized for the different algorithms in Section 5 including the calculation of Receiver Operating Characteristics (ROC),

• •

The detection sensitivity as a function of separation, referred to as contrast curve, is calculated in Section 6 where the contrast is defined as the companion to host star brightness ratio in a spectral band,

• •

The follow-up strategy and the vetting of point source candidates is discussed in Section 7,

• •

The contours of the planet completeness, which is the fraction of planets that could have been detected, are then derived in Section 8,

• •

We conclude in Section 9.

We refer any reader who is not familiar with the data processing of high contrast images to Appendix A, which includes a detailed description of KLIP, the matched filter and the planet PSF forward model. The mathematical notations are summarized in Appendix B.

2 Observations and Data Reduction

2.1 Observations

In this paper, we use 330 observations from the GPI Exoplanet Survey (GPIES) (Gemini programs GS-2014B-Q-500, GS-2015A-Q-500, GS-2015B-Q-500; PI: B. Macintosh) to construct and test our matched filter. A typical GPIES epoch consists in $38\times 1$ minute exposures in H band ($1.5-1.8\mu m$). The number of usable raw spectral cubes can be lower, due to degrading weather condition or isolated star tracking failures. We arbitrarily consider any dataset 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 bad conditions, which resulted in a number of datasets 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 datasets 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\times 281$ pixels in the spatial dimensions, half of which are however not filled with data due to the tilted IFS field of view. Therefore, a typical GPIES dataset includes approximately $1400$ images at different position angles and wavelengths.

2.2 Raw Data Reduction

Spectral data cubes are built from raw IFS detector images using standard recipes from the GPI Data Reduction Pipeline111Documentation available at http://docs.planetimager.org/pipeline/. version 1.3 and 1.4 (Perrin et al., 2016). The process includes correction for dark current, bad pixels, correlated read noise, and cryocooler vibration induced microphonics (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 solutions 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.

2.3 Use of Simulated Planets

Simulated planets are necessary to optimize and characterize a detection algorithm due to the scarcity of real point sources in high-contrast images. We decided to neglect for now the PSF smearing due to the sky rotation in a single exposure. The planets’ spectra are selected from atmospheric models described in Marley et al. (2017, in preparation) and Saumon et al. (2012). For a given cloud-coverage, the only model parameter having 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 spectrum is matching. 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 are cloudy objects whose spectra are dominated by H2O and CO. Of the known extrasolar planets, 51 Eridani b is an example of the T-type (Macintosh et al., 2015) and $\beta$ 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 spots spectrum with a stellar spectrum, which is then used to translate the spectrum from physical units to 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.

2.4 Speckle Subtraction with KLIP

Each individual image in the dataset is speckle subtracted using a Python implementation of KLIP called PyKLIP222Available under open-source license at https://bitbucket.org/pyKLIP/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 which can be selected from the same dataset 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-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 in small $100\ \text{pixel}$ arcs to which a $10\ \text{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 will be 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\arcsec$ ($\approx 71$ pixels).

3 Matched Filter

3.1 Concept

In the field of signal processing, a matched filter is the linear filter maximizing the SNR 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 dataset, each image is indexed by its exposure number $\tau$ and its wavelength $\lambda$. We define the vector ${\bm{p}}_{l}$, with $l=(\tau,\lambda)$, as a specific speckle subtracted image. Similarly, we define the matched filter template ${\bm{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 simplifies significantly 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 $\rho$ and position angle $\theta$ is then given by

[TABLE]

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 ${\bm{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 SNR of the planet can be written

[TABLE]

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

3.2 Matched Filter Template and Forward Model

The calculation of the matched filter template ${\bm{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 a 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 Section 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’s location on the image and therefore needs to be calculated for each pixel.

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 smaller values of the exclusion criterion. The test was performed with two representative datasets at two different separations: (a) HR 7012 at $0.3\arcsec$ and $0.6\arcsec$ and (b) HD 131435 at $0.4\arcsec$ and $0.8\arcsec$. The datasets 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 good dataset while HD 131435 is of median quality. A total of 35 simulated exoplanets were injected at each separation in seven copies of the same dataset. For each copy of the dataset, the five simulated planets, that are $72^{\circ}$ apart, were rotated by $10^{\circ}$ in position angle to better sample the image.

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 will be referred 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 dataset is collapsed before the matched filter is performed. It is outside 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).

3.3 Gaussian Cross Correlation

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 dataset is first derotated, coadded and then collapsed using the reduction spectrum. The resulting image is cross correlated with a 2-D Gaussian kernel with a Full-Width-Half-Maximum (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\ \text{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 SNR 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\times 20$ pixels stamp projected on a $281\times 281$ pixels image.

3.4 Gaussian Matched Filter

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 pixels 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\times 20$ pixels stamp from which a central disk with a $2.5\ \text{pixel}$ radius is removed. To summarize the Gaussian Matched Filter for GPIES, the template is a $20\times 20\times 37$ pixels cube stamp, which includes the spectral dimension, projected on a $281\times 281\times 37$ pixels datacube.

3.5 Forward Model Matched Filter

The Forward Model Matched Filter (FMMF) is performed on uncollapsed dataset 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-\Delta\rho,\rho+\Delta\rho]\times[\theta-\Delta\rho/\rho,\theta+\Delta\rho/\rho]$ and $\Delta\rho=10\ \text{px}$. In this case, the center of the arc is not masked, which results in an overestimation the standard deviation in the presence of a planet signal. This doesn’t 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 ${\bm{p}}_{l}^{\top}{\bm{m}}_{l}$, ${\bm{m}}_{l}^{\top}{\bm{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 SNR 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$ pixels multidimensional stamp, therefore including both spectral and time dimensions, projected on an uncombined $281\times 281\times 37\times 38$ pixels dataset.

The FMMF reduction of a typical GPI 38 exposure dataset requires around $30$ wall clock hours on a computer equipped with $32$ $2.3$ GHz cores and using $\sim 20$ GB of Random-Access Memory (RAM). As a consequence, a super computer is necessary to process an entire survey. Each dataset 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 of the order of $10^{6}$ CPU hours including the data processing necessary for the work presented in the remainder of the paper.

3.6 SNR Calculation

In practice, the theoretical SNR 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 SNR is overestimated due to overly optimistic assumptions on the noise distribution. The residual noise is mostly white and Gaussian in areas not dominated by speckle noise but both assumptions break down close to the mask due to speckle noise, causing the matched filter to lose some validity although it remains nonetheless relatively effective. 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 SNR. This is why an empirical standard deviation is estimated in the matched filter map as a function of separation using a $4$ pixels wide annulus as illustrated in Figure 5. In order to prevent a planet from biasing its own SNR, 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 a SNR relates to the calibrated SNR.

For a centered Gaussian noise, the SNR is related to the False Positive Rate (FPR) following

[TABLE]

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 a False Positive Rate equal to $\mathrm{FPR}=2.9\times 10^{-7}$, which represents a false detection every $3.4\times 10^{6}$ independent samples or equivalently of 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\,\mu m$.

3.7 Example Reduction

Figure 6 gives an example of reduction using the three different algorithms on the HR 7012 GPIES dataset observed on 2015 April 08. This dataset doesn’t contain any visible astrophysical signal. Simulated planets were injected according to Section 2.3. With the exception of the anomaly at $0.4\arcsec$, the FMMF consistently yields a higher SNR than the two other methods and the final image shows fewer residual features. In the following, the simulated datasets 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 it is not true with the Gaussian Matched Filter and the Gaussian Cross Correlation.

4 Optimization

The free parameters for KLIP and the matched filter are the exclusion criterion, the reduction spectrum, the number of Karhunen-Loéve modes, the number of reference images, and the shape of the sectors. In this paper, we are only optimizing with respect to the exclusion criterion and the spectrum of the template. The exclusion criterion is only optimized for T-type objects but we do not expect a significant difference for L-type objects. Preliminary tests show that the performance of the algorithm is less sensitive to the choice of the other parameters.

4.1 Exclusion Criterion

The exclusion criterion is defined in Section 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 SNR.

The optimal value of the exclusion criterion is found by calculating the SNR of simulated planets for different values of the parameter. We used similar simulated planets as for Figure 3 with a $1000\,\textrm{K}$ cloud-free model spectrum and a $600\,\textrm{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 optimimum around $0.7$. The optimal exclusion criterion doesn’t seem to significantly depend on separation or dataset quality. This apparent stability of the FMMF optimum suggests that a single value of the exclusion criterion can be used for the entire survey. Consequentially, 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 SNR than the two other algorithm but this does not necessarily mean it has better detection efficiency. Interestingly, neither the GCC nor GMF have a consistent or even a well-defined optimum. It is very common in the field to reduce a dataset 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 dataset.

4.2 Spectral Mismatch

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\,\mu m$ while L-type objects feature cloudy atmosphere with a flatter spectrum that peaks in the second half of the band.

We created ten copies of the HD 131435 GPIES dataset each with 16 simulated planets injected according to the spiral pattern from Figure 6. In each copy of the dataset, 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 does not significantly change the shape of the spectra. Then, each simulated dataset 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 SNR for the 16 simulated planets. Figure 9 shows the median SNR 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 weigh 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 outside the scope of this paper. We also conclude that only two spectra, the cloud-free $600\,\textrm{K}$ and the cloudy $1300\,\textrm{K}$, are necessary to allow the detection of most giant planets without a significant loss of SNR. This result is consistent with a similar study in Johnson-Groh (2017) using TLOCI (Marois et al., 2014). Although the $600\,\textrm{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\,\textrm{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\,\textrm{K}$ cloudy model spectrum as for the reduction.

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

5 Noise Characterization

5.1 Noise Distribution

By combining the $330$ GPIES H-band observations after removing any known astrophysical point sources, we have been able to estimate the Probability density Functions (PDF) of the residual noise up to an unprecedented precision in the fully reduced SNR 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 SNR 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 SNR in direct imaging. The Gaussian cross correlation has the same PDF as the SNR maps calculated from the speckle subtracted images with no cross correlation or matched filter. The Gaussian matched filter and the FMMF both significantly improve the statistics of the residuals but remain quite remote from an ideal Gaussian distribution. The excess of high SNR occurrences can be explained by either a truly non-Gaussian statistic or by a poor estimation of the standard deviation of the noise when calculating the SNR. For example, an underestimated standard deviation for a pixel will result in more high SNR false positives. The relative improvement between the cross correlation and the matched filter suggests that the latter explanation may still be dominant, which we will discuss in more detail in the next paragraph.

Figure 11 shows the spatial densities of $3\sigma$ false positives as a function of separation and position angle. Due to the correlation in the final SNR 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 SNR 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 SNR pixel is flagged, and a four pixel area radius is masked around it. This process is repeated until the only false positives left are below a predefined SNR threshold. Both the Gaussian Cross Correlation and the Gaussian Matched Filter 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 Gaussian Cross Correlation has a significant excess of false positives around $90^{\circ}$ and $270^{\circ}$ 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 North-East 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 Gaussian Matched Filter doesn’t 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 SNR detection threshold should yield the same number of false positives. The real performance of each algorithm will be 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 SNR threshold for each dataset and as a function of position but the lack of local independent samples to estimate a position dependent PDF at small false positive rates makes this endeavor very challenging. Indeed, one needs to have probed the PDF of the noise at high SNR 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 individual images.

The FMMF L-type reduction has significantly more false positives near the mask than a T-type reduction. This is likely due to the fact that 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 SNR.

5.2 Receiver Operating Characteristic

In general, an improvement in the SNR 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 SNR than the Gaussian matched filter 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; Jensen-Clem et al. 2017, submitted). 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. Multi-dimensional 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\,\textrm{K}$ cloud-free T-type spectrum, reduced with a $600\,\textrm{K}$ cloud-free model spectrum, or a $1300\,\textrm{K}$ cloudy L-type spectrum, reduced with a the same cloudy spectrum. All planets were injected at a $4\times 10^{-6}$ contrast in 330 GPIES H-band datasets using the spiral pattern illustrated in Figure 6. The advantage of using a large number of datasets 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. Also, the SNR threshold corresponding to a given false positive rate is always higher in the case of the cross correlation than for the matched filters due to the larger tail in the PDF. For example, it has a T-type false positive rate roughly $6$ times higher at $5\sigma$ than the two other algorithms.

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 SNR threshold therefore depends on the algorithm as shown in Figure 12.

6 Contrast Curve

In this paper, contrast is defined as the broadband flux ratio between the companion and star. The contrast curve is defined as the $50\%$ detection completeness contour assuming a false positive rate sufficiently low to limit the number of false positives. The false positive rate can be expressed in terms of a SNR threshold, which is not necessarily $5\sigma$. Using a hard $5\sigma$ threshold does not allow for a meaningful comparison of contrast curves, because different algorithms can lead to different numbers of false positives for the same SNR (as demonstrated in Section 5.2). Contrast curve calculations require a calibration step to translate pixel values to planet contrast, which is in some cases referred to as throughput correction. Indeed, speckle subtraction algorithm like KLIP are known to over-subtract the signal of the planet and make it appear fainter than it really is. This effect can be calibrated out by inflating the standard deviations by a certain factor, known as throughput, calculated from simulated planet injection. However, the throughput correction only makes sense if the pixel values, and therefore the standard deviation, are in units of contrast, which is not the case in this paper. The matched filter maps don’t directly estimate the contrast of a planet but rather try to estimate a theoretical SNR as defined in Equation 2. Therefore, the well known throughput correction is here replaced by a conversion factor to relate the matched filter values to actual contrast estimates, which is also calibrated using simulated planets.

As a consequence, we define a contrast curve as $\eta\gamma(\rho)\sigma_{\mathcal{S}}(\rho)$ where $\eta$ is the SNR threshold, $\gamma(\rho)$ is a median conversion factor between the matched filter map and the true contrast of the planet and $\sigma_{\mathcal{S}}(\rho)$ is the standard deviation of the noise in the matched filter map.

The different thresholds $\eta$ were determined in Section 5.2 for the three algorithms to yield a false positive frequency of $0.05$ per epoch. The standard deviation $\sigma_{\mathcal{S}}(\rho)$ is calculated in concentric annuli according to Section 3.6. The conversion factor $\gamma$ is empirically determined by injecting simulated planets with known contrast in each dataset. The contrast of the simulated planets is chosen to result in a signal somewhere between $5-15\sigma(\rho)$ in the final matched filter map. In total, 128 planets are injected at 16 different separations and position angles in eight different copies of the dataset. The original matched filter map is subtracted from the simulated planet reductions in order to remove the effect of the residual noise. The conversion factor is linearly interpolated from the median of the eight simulated planets at each separation as shown in Figure 13. Note that the contrast curve is a function of the spectrum of the simulated planets as well as of the reduction spectrum. Ideally, there should be as many contrast curves as there are possible spectra. We have limited our study to objects with either a 1000K cloud-free T-type spectrum, reduced with a 600K cloud-free model spectrum, or a 1300K cloudy L-type spectrum, reduced with the same cloudy spectrum. A caveat is that the contrast curves are therefore only truly valid for planets with the same spectrum as the injected planets.

The GPIES median contrast curve for each algorithm as well as the associated detection threshold for both spectral type reduction is given in Figure 14. FMMF yields the best median contrast curve at all separations. In the T-type reduction and compared to the Gaussian Matched Filter, the FMMF contrast enhancement ranges from a median $25\%$ up to more than a factor 2 in some cases. The L-type median contrast enhancement drops from $25\%$ below $0.5\arcsec$ to $10\%$ at larger separation. The median contrast improvement relative to the cross-correlation is around $50\%$, which corresponds to a factor 2.3 gain in exposure time assuming a square root increase of the SNR with time. Tests have shown that this assumption generally holds true to the exception of a slight dependence to observing conditions. The local maximum in the contrast curve for the Gaussian Cross correlation and the Gaussian matched filter is likely due to the non-optimal definition of the reduction sectors for a regular KLIP reduction. We will therefore ignore the FMMF contrast increase in the $0.3\arcsec-0.4\arcsec$ range as we do not consider it representative of the overall performance.

The sensitivity of any observation is highly dependent on the observing conditions and the brightness of the star. Figure 15 shows the percentiles contours of the contrast curves for GPIES using FMMF. Figure 16 shows the FMMF median contrast curves as a function of the I magnitude of the star. There is an order of magnitude ratio between the sensitivity around the faintest stars ($I_{\textbf{mag}}\approx 9$) of the survey compared to the brightest stars ($I_{\textbf{mag}}\approx 2$).

7 Candidate vetting

In order for the contrast curves and the planet completeness to be accurate, any signal above the detection threshold must be properly vetted. This means that any high SNR signal should be confirmed as a true or a false positive. The best approach is always to follow up all the candidates, which can be a significant telescope time investment. In the case of a real astrophysical signal, the second epoch is generally used to determine proper motion and parallax in order to exclude or confirm the possibility of a background object. In the case of low SNR candidates, if it cannot be detected in the second epoch, it is usually classified as a false positive. However, it is then necessary to improve the contrast curve in the second epoch in order to exclude the real signal hypothesis with a high enough significance. We discuss the necessary contrast improvements in this section.

First, the detection threshold can be lowered in the follow-up observation compared to the first epoch, because the constraint on the noise is higher. A spurious signal would need to be in the same position in order for it to be mistakenly considered as a true detection. With a $4\sigma$ detection threshold in the second epoch, there is less than a chance in a thousand to wrongly classify the first detection as a real signal, which seems reasonable. Indeed, there is on the order of $10^{3}$ independent samples in a $1\arcsec$ field of view in a GPI image assuming a conservative characteristic correlation length of the residuals of 3 pixels and a $4\sigma$ threshold yields less than a false positive per epoch for both matched filters.

Second, we can estimate the probability of detecting the signal in the new epoch as a function of the original SNR and the ratio of contrast standard deviation $\sigma_{\epsilon 1}/\sigma_{\epsilon 2}$. Assuming a known contrast for a point source, the probability of detection is given by the tail distribution of the planet signal at the detection threshold in the second epoch. However, it needs to be marginalized over the planet contrast due to its uncertainty in the first epoch. The detection probability in the second epoch assuming a detection threshold $\eta_{2}$ can be written as

[TABLE]

with $\mathrm{PDF}$ and $\mathrm{CDF}$ respectively the Probability Density Function and the Cumulative Distribution Function of a Gaussian distribution with zero mean and unit standard deviation. Figure 17 shows the detection probability contours as a function of the original SNR and the noise ratio.

The problem is that candidate follow-ups are expensive in telescope time and a significant fraction of the detections are expected to be false positives. Indeed, a significant number of candidates appear to be reduction artifacts. Some of these artifacts are very sensitive to the edge of the sectors and can therefore be spotted by running a new reduction with an annulus centered at their separation. Real objects are less sensitive to the definition of the sectors so if the SNR of the candidate drops significantly, it is likely a false positive and should not be reobserved. Figure 18 shows the L-type candidates above $5.2\sigma$ as a function of their original SNR and the relative SNR with the new reduction. Because the SNR ratio for the real point sources are always close to unity, we conclude that we can reject all candidates for which this ratio drops under 0.8 without significantly biasing the algorithm completeness. This plot cannot be shown for the T-type reduction because there is only one confirmed T-type object in GPIES, which is 51 Eridani b. Most point sources are indeed background stars and the few substellar objects feature a L-type spectrum. The lack of control sample for the T-type reduction does not allow us to define a boundary in this case. We acknowledge that rejecting candidates based on an ad-hoc criterion is not the best solution and proves that work remains to be done. Ideally, we would want to combine the results of both reductions into a single detection metric rather than applying different cutoffs sequentially.

8 From Contrast Curve to Completeness

The main goal of a direct imaging exoplanet survey, after discovering new worlds to characterize, is to put constraints on the population of wide orbit sub-stellar objects (Nielsen et al., 2013; Bowler et al., 2015; Galicher et al., 2016; Vigan et al., 2017). Deriving the frequency of exoplanets first relies on the completeness of the survey. The most common method to calculate the completeness of an observation is to reduce the same dataset with simulated planets at all separations and contrasts. The completeness at a given point is then the fraction of the fiducial planets detected. This approach requires many reductions for each dataset, which is tractable with classical ADI reductions but computationally out of reach for the FMMF. Lafrenière et al. (2007b) showed that, assuming a Gaussian noise, the completeness can be estimated directly from the contrast curve. For example, a planet lying exactly on the contrast curve will have a $50\%$ chance to fall above or below it due to the noise. The probability of detecting a planet of a given contrast is given by the tail distribution of a Gaussian distribution centered on the contrast of the planet evaluated at the contrast curve. We have shown that the Gaussian assumption is not verified in the tail of the noise distribution but this effect would only become important for completeness of planets far from the contrast curve. In the ideal case, the $(\eta+1)\gamma(\rho)\sigma_{\mathcal{S}}(\rho)$ and the $(\eta-1)\gamma(\rho)\sigma_{\mathcal{S}}(\rho)$ curves respectively represents the $84\%$ and $16\%$ completeness contour. One caveat is that the azimuthal variations of the conversion factor (e.g. due to the wind-butterfly) add a scatter term to the measured contrast of the planet effectively widening the completeness contours. For example, a planet inside the wind-butterfly suffers from a higher conversion factor than a planet outside of it, resulting in a lower detection probability.

The azimuthal variations of the conversion factor due to the wind-butterfly or other artifacts can be interpreted as a noise term with a standard deviation $\sigma_{\gamma}$. Assuming that the residual noise and the conversion factor variations are independent, the standard deviation of the contrast can be written

[TABLE]

which is used to estimate the completeness. Figure 19 shows two examples of completeness contours, one with a strong wind-butterfly and another where it is negligible.

9 Conclusion

In order to improve GPIES planet sensitivity, we implemented a new matched filter based algorithm using KLIP and a forward modeled PSF template. This algorithm includes the speckle subtraction efficiency of KLIP, while mitigating the PSF distortion penalty. We also presented a Gaussian cross correlation and a Gaussian matched filter for comparison. The cross correlation requires an identically distributed noise in order to be valid. This assumption is not verified in GPI images and it therefore leads to an increase of false positives in the part of the image, like the wind-butterfly, where the noise is consistently stronger. An accurate matched filter needs to include a normalization by the local noise standard deviation, which improves significantly the statistics of the residual noise in the final SNR maps.

After optimizing the aggressiveness of the algorithm, we showed that two reduction spectra, a T-type and a L-type, are sufficient to recover sub-stellar companions candidates with their expected spectra without a significant loss in SNR.

The uniform reduction of more 330 datasets from GPIES allowed the unprecedented characterization of the PDF of the residual noise, which significantly deviates from a Gaussian distribution. We showed that the FMMF has a uniform spatial distribution of false positives that does not feature any excesses due to the sector boundaries or the wind-butterfly. The improved performance provided by the FMMF also appeared in the ROC curves and the contrast curves that have been shown for GPIES. We built the contrast curves using a different detection threshold for each spectral type and algorithm such that it yields the same number of false positives in each case. The FMMF allows for example the detection of objects that are $50\%$ fainter than a classical cross correlation. We reaffirm that a hard $5\sigma$ detection threshold doesn’t allow meaningful comparison of algorithms performance and can overestimate the sensitivity of a survey. In the future, more diverse algorithms should be compared following a similar recipe to the one presented in this paper.

To finish, we showed how the planet completeness can be derived from each contrast curve, which is simply the $50\%$ detection completeness contour. These completeness contours can be used in a future exoplanet population study.

Acknowledgements: This work is based on observations obtained at the Gemini Observatory, which is operated by the Association of Universities for Research in Astronomy under a cooperative agreement with NSF on behalf of the Gemini partnership, whose membership includes: NSF (United States); the National Research Council (Canada); the Comisión Nacional de Investigación Científica y Tecnológica (Chile); the Australian Research Council (Australia); the Ministério da Ciência, Tecnologia e Inovação (Brazil); and Ministerio de Ciencia, Tecnología e Innovación Productiva (Argentina). GPI data are archived at the Gemini Observatory Archive333http://www.gemini.edu/sciops/data-and-results/gemini-observatory-archive : https://archive.gemini.edu/searchform. The research was supported by grants from NSF, including AST-1411868 (J.-B.R. B.M., K.F., A.R., J.L.P.) and AST-1518332 (J.J.W., R.J.D, J.R.G., and P.K.). Support was provided by grants from NASA, including NNX14AJ80G (B.M., E.N., and V.P.B.), NNX15AD95G (J.J.W., R.J.D, J.R.G., and P.K.). This work benefited from NASA’s Nexus for Exoplanet System Science (NExSS) research coordination network sponsored by NASA’s Science Mission Directorate. K.M.M.’s and T.B.’s work is supported by the NASA Exoplanets Research Program (XRP) by cooperative agreement NNX16AD44G. M.A.M.B. was supported for this work by NASA through Hubble Fellowship grant #51378.01-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. J.R and R.D acknowledge support from the Fonds de Recherche du Quebec.

Appendix A Method

In this appendix, we summarize the key theoretical results used in this paper. First, the Karhunen-Loève Image Processing (KLIP) framework and its effect on the planet PSF, i.e. self-subtraction, is introduced in Appendix A.1. Then, the principles of a matched filter are presented in Appendix A.2. Finally, the forward model is defined in Appendix A.3.

Throughout the paper, with some rare exceptions, scalars are lower case Greek characters, vectors are bold lower case Latin characters and matrices are bold upper case Latin characters. A summary of the notations can be found in Appendix B.

A.1 Karhunen-Loève Image Processing

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 filter 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 will 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 Section A.1.2. A detailed description of the mathematical formalism underlying the KLIP algorithm can be found in Savransky (2015), which we quickly 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

[TABLE]

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

[TABLE]

with its eigenvectors and eigenvalues ${\bm{v}}_{1},{\bm{v}}_{2},\dots,{\bm{v}}_{N_{R}}$ and $\mu_{1},\mu_{2},\dots,\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 square sense. They are defined as

[TABLE]

The images can be written in a more compact form by defining the matrices $\boldsymbol{V}_{K}=[{\bm{v}}_{1},{\bm{v}}_{2},\dots,{\bm{v}}_{K}]^{\top}$, $\boldsymbol{Z}_{K}=[{\bm{z}}_{1},{\bm{z}}_{2},\dots,{\bm{z}}_{K}]^{\top}$ and the diagonal matrix $\boldsymbol{\Lambda}=\textrm{diag}\left(\mu_{1},\mu_{2},\dots,\mu_{K}\right)^{\top}$, with $K\leq 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}=[{\bm{z}}_{1},{\bm{z}}_{2},\dots,{\bm{z}}_{K}]^{\top}$ is equivalently written as

[TABLE]

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

[TABLE]

where $p$ refers to the processed image. The science and the reference images are often part of the same dataset and the science image simply 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, however, 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 A5, following Pueyo (2016). In this section only, variables defined in Section A.1.1 are assumed to be planet free, and the hat indicates the perturbed variable, such that

[TABLE]

with $\Delta{\bm{r}}_{j}=\epsilon{\bm{a}}_{j}$, where ${\bm{a}}_{j}$ is the normalized planet signal. The planet signal is defined by its amplitude $\epsilon$, 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 $\epsilon$ is in practice the contrast of the planet. Similarly, we write $\hat{{\bm{z}}}_{j}={\bm{z}}_{j}+\Delta{\bm{z}}_{j}$, $\hat{\boldsymbol{Z}}_{K}=\boldsymbol{Z}_{K}+\Delta\boldsymbol{Z}_{K}$, $\hat{{\bm{i}}}={\bm{i}}+\Delta{\bm{i}}={\bm{i}}+\epsilon{\bm{a}}$, and $\hat{{\bm{p}}}={\bm{p}}+\Delta{\bm{p}}\approx{\bm{p}}+\epsilon{\bm{m}}$, where ${\bm{m}}$ is the planet model that will be used to build the matched filter template ${\bm{t}}$ introduced in Equation 2.

Applying KLIP to the perturbed images gives

[TABLE]

and expanding LABEL:eq:klipPert, one identifies ${\bm{p}}$, which should only contain residual speckles such that

[TABLE]

Then, the perturbed processed image is given by

[TABLE]

The second and third terms of LABEL:eq:selfsubtraction are the first order distortions introduced by KLIP on the planet signal as shown in Figure 2 (b). The over-subtraction 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 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 $\Delta\boldsymbol{Z}_{K}=0$, which is the case when using a Reference star Differential Imaging (RDI) approach. Using the Locally Optimized Combination of Images (LOCI) Lafrenière et al. (2007a), 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 possible 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 overlap in the reference images and in the science image spatially. The distortion of the central peak of the PSF is more important when the overlap is big. 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 a the relative flux overlap with the planet signal in 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 accepts more or less images in the reference library when the wavelengths are different. This implementation of the exclusion criterion has given significantly better result than any of the two other methods. The conversion between the three metrics is given for a specific planet PSF flux ratios 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 dataset. 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.

The choice of the exclusion criterion is a trade-off between an effective speckle subtraction and a limited self-subtraction. A small value of the exclusion criterion includes more correlated images in the reference library and results in a more aggressive reduction, due to the larger degree of self-subtraction. The aggressiveness is also affected by the number of Karhunen-Loève modes used in the subtraction, however 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.

A.2 Matched Filter

In the field of signal processing, a matched filter is the linear filter maximizing the SNR 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 dataset 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_{\lambda}$ is the number of spectral channels, and $N_{\tau}$ is the number of exposures. A typical GPIES dataset is the concatenation of all the processed images $\left\{p_{l}\right\}$ 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 dataset contains a single point source, it can be decomposed as

[TABLE]

With ${\bm{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 $\epsilon$ is again the contrast. It is the concatenation of the planet model PSF for each image $\left\{{\bm{m}}_{l}\right\}$ 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 ${\bm{\Sigma}}$. For conciseness, the $(\rho,\theta)$ notation in ${\bm{t}}(\rho,\theta)$ will be dropped in the following 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:

[TABLE]

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

[TABLE]

The theoretical SNR of a point source is given as a function of position by,

[TABLE]

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.

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

[TABLE]

where $\sigma_{l}$ is the local standard deviation at the position $(\rho,\theta)$. Note that the planet model ${\bm{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.

The theoretical SNR of the planet becomes

[TABLE]

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

A.3 KLIP Forward Model

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 due to 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 $\frac{\Delta p}{\epsilon}$ from LABEL:eq:selfsubtraction is what should be used as the template for the matched filter (a.k.a. planet signature in Cantalloube et al. (2015) and 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 $\frac{\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 $\Delta\boldsymbol{Z}_{K}$ because it is a nonlinear function of the planet flux. This is why 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:

[TABLE]

With $\boldsymbol{A}=\left[{\bm{a}}_{1},{\bm{a}}_{2},\dots,{\bm{a}}_{N_{R}}\right]^{\top}$ and $\boldsymbol{C}_{\boldsymbol{A}\boldsymbol{R}}=\boldsymbol{A}\boldsymbol{R}^{\top}+\left(\boldsymbol{A}\boldsymbol{R}^{\top}\right)^{\top}$.

The forward model of the planet, a.k.a template, is then given by

[TABLE]

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

Appendix B Mathematical Notations

We present here a summary of the notations used in this paper. Table 1, Table 2, and Table 3 respectively contain the definition of the scalars, the vectors and the matrices.

\startlongtable