© O. Flasseur et al. 2020

1. Introduction

Direct imaging from the Earth is a method of choice for the detection and characterization of exoplanets (Traub & Oppenheimer 2010). By confronting numerical models of estimated astrometry and photometry (Vigan et al. 2010), direct imaging gives access to quantitative physical properties on the detected sources such as age, mass, and effective temperature (Chabrier et al. 2000; Allard et al. 2003, 2007). The most efficient observation strategy is based on angular differential imaging (ADI, see Marois et al. 2006) which consists of tracking the observation target over time (the telescope pupil remaining stable while the whole field of view rotates). A temporal sequence of images is thus acquired in which the companions describe a predictable motion in time while the speckle background remains quasi-static. This observation method is used routinely on the main instruments dedicated to the direct detection of exoplanets such as VLT/SPHERE (Beuzit et al. 2019; Chauvin et al. 2017), GEMINI/GPI (Macintosh et al. 2014; Greenbaum et al. 2018), Keck/VORTEX (Howard et al. 2010; Mawet et al. 2019), Magellan/MagAO (Morzinski et al. 2014; Wagner et al. 2018), and Subaru/SCExAO (Jovanovic et al. 2015; Currie et al. 2017). The recorded ADI images are then combined by a data-processing algorithm to disentangle the signature of the exoplanets from the speckles (Lafrenière et al. 2007; Marois et al. 2013, 2014; Soummer et al. 2012; Smith et al. 2009; Cantalloube 2016; Gonzalez et al. 2016, 2018; Flasseur et al. 2018a). The statistical modeling of the data plays a central role in these processing algorithms to control the probability of false alarms and prevent spurious detections (Mawet et al. 2014).

The main constraint of this method is due to the very high contrast between the target star and the exoplanets (≥10⁵ in infrared). In addition, the detection performance and the achievable contrast are strongly dependent on both the total parallactic rotation of the field of view and the quality of the observations often impacted by several artifacts. These artifacts can be spatially localized (e.g., in case of defective pixels) or can impact a larger fraction of the field of view in the form of large fluctuations when a decentering of the coronagraph or a sudden degradation of the adaptive optics (AO) correction occurs (e.g., a low-wind effect, Sauvage et al. 2015; Milli et al. 2018). The latter two effects are especially problematic since the exoplanets are most likely to be located close to the center of the field of view which is the area most affected by the fluctuations.

To the best of our knowledge, there is no particular strategy implemented in the state-of-the-art post-processing algorithms to deal with the evolution of the local quality of the acquired frames. Images presenting large fluctuations compared to the others are simply flagged as “bad frames” and discarded from the ADI stack, even if some areas of these images contain useful information. In this paper, we present an extension of our recently proposed detection algorithm accounting for the patch covariances of the background (PACO algorithm, see Flasseur et al. 2018a). The proposed improvement is based on the modeling of the temporal variations of the amplitude of the background fluctuations jointly to the spatial covariances. To improve the robustness of the method, we spatially weight each temporal frame in a data-driven fashion based on its relative degree of fluctuation.

Section 2 describes our local modeling of the spatio-temporal fluctuations and Sect. 3 details the extension of the PACO algorithm. In Sect. 4, we use VLT/SPHERE-IRDIS datasets to illustrate the performance of the proposed method in terms of detection capability, achievable contrast, and astrometry accuracy in comparison to two state-of-the-art algorithms routinely used on ADI datasets, TLOCI (Marois et al. 2013, 2014) and KLIP (Soummer et al. 2012), both implemented in the SpeCal software (Galicher et al. 2018) of the SPHERE data center (Delorme et al. 2017). We show that by preventing the suppression of bad frames, the automatic local weighting of the images improves the achievable contrast on the whole field of view.

2. Local modeling of spatio-temporal fluctuations

Even in the absence of sources in the field of view other than the star masked out by the coronagraph, there are fluctuations in images from an ADI sequence. There are several factors that induce these fluctuations, in particular: (i) evolution of the point spread function (PSF) due to changes in the observing conditions and alterations of the quality of the adaptive optics correction; (ii) photon and thermal noise on the camera as well as the erratic response of some uncorrected bad pixels; (iii) partial decentering of the coronagraph; (iv) and evolving noncommon-path aberrations (uncorrected by the AO).

In Fig. 1a, we illustrate the temporal evolution of the intensity measured with the SPHERE-IRDIS instrument. We selected some frames from an observation of HIP 72192 (out of T = 96 temporal frames acquired under the 2015-06-11 – 095.C-0389 ESO program with a total apparent rotation of the field of view of 17.3°) that displayed particularly strong fluctuations due to difficult observing conditions. Since these fluctuations are mostly located around the coronagraph, we show a zoom on the central region of the images corresponding to a field of 2.8 × 2.8 arcsec². The fluctuations are spatially structured and their magnitude varies from one image to the other in the ADI sequence.

Fig. 1. Central region of SPHERE-IRDIS images of HIP 72192: panel a: measured intensity for five selected frames showing strong temporal fluctuations; panel b: map of the estimated temporal scaling parameters ${\widehat{\sigma }}_{k,\ell }$ for the matching frames: ℓ = 1, 3, 8, 13 and 14.

In Flasseur et al. (2018a), we modeled the statistical distribution of the patch r_θk, tℓ: the small window made of a few tens of pixels centered at angular location θ_k in the ℓ-th frame. This patch with K pixels was described as the possible superimposition of the signal of point-source (if present) and a background patch f_θk, tℓ, considered as a random realization of a multivariate Gaussian 𝒩(m_k, C_k), with spatially variant mean m_k ∈ ℝ^K and covariance C_k ∈ ℝ^K × K. At any given location θ_k, the background fluctuations captured through the covariance matrix C_k were considered stationary (i.e., not evolving with time).

To account for the presence of “bad” frames with stronger background fluctuations than the other frames, we consider in this paper that the background patch f_θk, tℓ is a random realization of a Gaussian scale mixture (GSM; see e.g. Wainwright & Simoncelli 2000), also known as a compound-Gaussian model (see e.g. Conte et al. 1995):

$$\begin{array}{c}\hfill {\mathit{f}}_{{\theta }_k,t_{\ell }}={\mathit{m}}_k+{\sigma }_{k,\ell }{\mathit{u}}_{k,\ell },\end{array}$$(1)

where σ_k, ℓ >  0 is a scalar (random or deterministic) variable that acts as a scaling parameter and the vectors u_k, ℓ ∼ 𝒩(0, C_k) are independent centered Gaussian random vectors.

Gaussian scale mixtures cover a wide range of distributions; in particular heavy-tailed distributions like gamma distribution, Student’s distribution, or the generalized Laplacian. They can account for the presence of outliers in the data, such as large deviations. They have been introduced in image processing to model the distribution of wavelet coefficients (see Wainwright & Simoncelli 2000; Portilla et al. 2003) and are widely used to model heterogeneous clutter in synthetic aperture radar; see for example Conte et al. (1995), Gini & Greco (2002), Wang et al. (2006), Pascal et al. (2008).

Rather than directly handling the GSM distribution, we estimate the scale parameter σ_k, ℓ for the patch centered on θ_k that is extracted from frame t_ℓ, and consider that f_θk, tℓ is distributed according to $\mathcal{N}({\mathit{m}}_k,{\sigma }_{k,\mathcal{l}}^2{\mathit{C}}_k)$. This corresponds to considering parameter σ_k, ℓ as a deterministic nuisance parameter; see Pascal et al. (2008).

Figure 1b displays maps of the estimated scaling factors ${\widehat{\sigma }}_{k,\ell }^2$ for the five frames shown in Fig. 1a. In these maps, large values of the scaling factors match the areas in which the intensities are much larger than the average. The variations observed within the field of view indicate that temporal fluctuations are spatially nonstationary: abnormally large fluctuations occur only in some areas and cannot be compensated for by a factor common to the whole image. By closer visual inspection, the scaling factors can be seen to vary over very small distances. This indicates why methods based on linear combinations of template on-axis PSFs (TLOCI) or modes (KLIP) fail to capture such small-scale variations, even when a local fitting based on annuli or angular sectors is performed.

An analysis of the empirical distribution of intensities in a patch, in the absence of a source, is carried out in Fig. 2. Three cases are compared: (i) patches extracted at a location θ_k close to the coronagraph, selected to highlight temporal fluctuations (first two rows of the figure); (ii) patches extracted at a location θ_k farther from the coronagraph, showing only moderate temporal fluctuations (following two rows); and (iii) patches extracted all over the field of view (last row). The temporal collection of patches {r_θk, tℓ}_{ℓ = 1 : T} is considered in each case. The empirical distribution of the intensities of the centered patches ${\{{\mathit{r}}_{{\theta }_k,t_{\ell }}-{\widehat{\mathit{m}}}_k\}}_{\ell =1:T}$ is displayed in the left column. Based on a local estimate ${\widehat{\mathbf{C}}}_k$ of the spatial covariance matrix, the distribution of spatially whitened and centered patches is then displayed on the central column. The Cholesky factor¹ ${\widehat{\mathbf{L}}}_k$ (i.e., the lower-triangular matrix such that ${\widehat{\mathbf{L}}}_k{\widehat{\mathbf{L}}}_k^{\mathrm{t}}={\widehat{\mathbf{C}}}_k^{-1}$) is used to obtain the collection of spatially whitened patches ${\{{\widehat{\mathbf{L}}}_k^{\mathrm{t}}({\mathit{r}}_{{\theta }_k,t_{\ell }}-{\widehat{\mathit{m}}}_k)\}}_{\ell =1:T}$. The empirical distribution of that collection is plotted. The right column displays the distribution of patch intensities when both a spatial whitening and a time-specific scaling is performed. The empirical distribution of the collection ${\{\frac{1}{{\widehat{\sigma }}_{k,\ell }}{\widehat{\mathbf{L}}}_k^{\mathrm{t}}({\mathit{r}}_{{\theta }_k,t_{\ell }}-{\widehat{\mathit{m}}}_k)\}}_{\ell =1:T}$ is plotted in this latter case.

Simply removing the temporal mean does not lead to independent and identically distributed normal residuals, as can be observed in the left column of Fig. 2. The variance of the residuals is much larger near the coronagraph than farther away, a behavior that is generally addressed by introducing a radial scaling of the residuals (Mawet et al. 2014). Locally modeling the spatial correlations provides a good approximation of the distribution of the patches in most locations θ_k, but fails in regions with large temporal heterogeneity: the distribution in that case is then better described by a GSM; see the central column of Fig. 2b. The presence of frames with larger fluctuations leads to an over-estimation of the variances in ${\widehat{\mathbf{C}}}_k$, which in turn leads to whitened samples with a variance that is less than one (Fig. 2e). When both a local modeling of the spatial covariance and a time-specific local scaling factor are considered, the empirical distributions are well modeled: after whitening, they closely match a standard Gaussian in all cases (near the coronagraph, farther from the coronagraph, and on average over the whole field of view). Our modeling and the local patch-whitening transform that derives from the GSM model can be seen as an extension of the radial-scaling strategy that is refined in several ways: (i) locality (vs. identical processing of an annulus), (ii) modeling of spatial correlations (vs. uncorrelated noise assumption), and (iii) modeling of time fluctuations (vs. constant correction).

Fig. 2. Temporal evolution and empirical distribution of the intensities within a patch: (a)–(b) in patches extracted at a location θk close to the coronagraph; (c)–(d) in patches extracted at a location θk farther from the coronagraph; (e) in all the patches from the field of view. Rows a and c: values within the patches for ten specific frames: t1, t3, t8, t13, t14, t25, t40, t54, t63 and t69. Rows b, d, and e: empirical histograms computed over all frames t1 to t69 and, in dashed line, a standard Gaussian. The first column corresponds to centered patches ${\{{\mathit{r}}_{{\theta }_k,t_{\ell }}-{\widehat{\mathit{m}}}_k\}}_{\ell =1:T}$. The second column corresponds to centered patches, after whitening the spatial correlations. The last column corresponds to centered patches that have been both whitened for the spatial correlations and equalized with the temporal scaling factors.

3. Adaptation of the PACO algorithm

3.1. Estimation of the statistics of the background

Estimating, at each spatial location, the mean background m, the spatial covariance C, and the scale parameters σ₁ to σ_T is very close to the problem of covariance structure estimation in radar under a compound-Gaussian clutter model (Pascal et al. 2008). There are two differences: (i) in radar the signal and covariance matrices are complex-valued; and (ii) the mean value is zero.

We derive in Appendix A the expression of the maximum likelihood estimates for the scale parameters, the local mean and the local covariance. The expression of each estimator depends on the others and so a fixed point is sought to obtain the joint estimator of the local background statistics. Based on the alternating application of formulas (A.2)–(A.4), an iterative algorithm is obtained. The algorithm is detailed in boxes 1 and 2.

Up to the centering, on accounting for the nonzero mean, the obtained fixed-point algorithm matches the estimator derived in Conte et al. (2002) in the context of clutter with a deterministic texture, and also the approximate maximum likelihood estimator derived in Gini & Greco (2002) under a stochastic model of texture. The convergence of the fixed-point algorithm has been established in Pascal et al. (2008) in the real-valued case as well as in the complex-valued case.

Since our estimates are based on small sample sizes, we include a shrinkage procedure similar to that proposed in Chen et al. (2010) and Flasseur et al. (2018a). To account for the weighting of the samples, we compute an effective number of frames $\stackrel{\sim }{T}={(\sum _{\ell }1/{\widehat{\sigma }}_{\ell }^2)}^2/(\sum _{\ell }1/{\widehat{\sigma }}_{\ell }^4)$ that corresponds to the variance reduction reached when performing weighted means with the weights $1/{\widehat{\sigma }}_{\ell }^2$. This effective number of frames is used to compute the shrinkage factor in Algorithm 3. This shrinkage operation can either be applied at each iteration (during step 4 of Algorithm 1) or after the alternating updates have converged. We tested both approaches and found that they converged to solutions leading to similar detection performances. Applying the shrinkage at each iteration ensures that matrix $\widehat{\mathbf{C}}$ never becomes singular or ill-conditioned.

Figure 3 illustrates the fast convergence of the alternating scheme. Each graphic corresponds to a different location θ_k in the field of view, depicted by a red dot. From left to right, they correspond respectively to a region with good temporal stationarity, medium temporal stationarity, and low temporal stationarity (close to the coronagraph). The convergence of Algorithm 1 is assessed based on 1000 random initializations of the weights ${\widehat{w}}_{\ell }$ (following a uniform distribution on [0,1]). The vector of weights at convergence when starting from a constant vector of weights ($\forall \ell ,{\widehat{w}}_{\ell }=1/T$) is used as a reference. At each iteration, the Euclidean distance to the reference vector of weights is reported, with normalization by the largest distance of all the random draws. Very fast convergence is observed to a unique solution, in all random trials. This empirical evidence indicates that the algorithm also converges when the shrinkage step is included at each iteration.

Fig. 3. Illustration of the convergence of Algorithm 1: 1000 random initializations of the weights ${\{{\widehat{w}}_{\ell }\}}_{\ell =1:T}$ converge to the same solution as that obtained by initializing with constant weights. Convergence at three locations θk depicted on the field of view is shown, covering three different cases: small dispersion of the value of scaling parameters (left), medium dispersion (center), and large dispersion close to the coronagraph (right).

3.2. Robust computation of a detection map

The PACO algorithm, introduced in Flasseur et al. (2018a), produces a detection map based on a hypothesis test evaluated at each pixel of the field of view. The ADI sequence of images is processed at the scale of image patches. The size of the patches is selected to capture the core of the off-axis PSF, for example a disk of K = 49 pixels for SPHERE-IRDIS in K1-K2 observing mode. If a point source object is located at angular position ϕ₀ in some reference frame, then its angular location ϕ_ℓ in the image at time t_ℓ can be deduced from the telescope pointing information by accounting for the rotation of the whole field of view during ADI acquisitions. Source detection is based on the analysis of the collection {r_{⌊ϕℓ⌉, tℓ}}_{ℓ = 1 : T} of the T patches formed by tracking the source through its apparent motion: each patch is extracted at the expected location, rounded to the closest pixel ⌊ϕ_ℓ⌉, of the point source with reference coordinate ϕ₀.

During the detection phase, we assume that the observations are dominated by the background signal due to stellar leakages and by the noise so that background statistics can be computed directly from the data (off-axis point sources are considered negligible at that step). The local statistics of the background are then computed with Algorithm 1, for each location ⌊ϕ_ℓ⌉, based on the collections {r_{⌊ϕℓ⌉, t}}_{t = 1 : T} of observed patches all centered at that same location ⌊ϕ_ℓ⌉. Once the mean ${\widehat{\mathit{m}}}_{\lfloor {\varphi }_{\ell }\rceil }$, spatial covariance ${\widehat{\mathbf{C}}}_{\lfloor {\varphi }_{\ell }\rceil }$ and scaling factor² ${\widehat{\sigma }}_{\lfloor {\varphi }_{\ell }\rceil ,\ell }$ have been estimated for all locations ⌊ϕ_ℓ⌉ corresponding to the trajectory of a hypothetical point source with reference location ϕ₀, the likelihood of two hypotheses can be compared:

$$\begin{array}{c}\hfill \{\begin{array}{cc}{\mathcal{H}}_0:{\{{\mathit{r}}_{\lfloor {\varphi }_{\ell }\rceil ,\ell }\}}_{\ell =1:T}\hfill & ={\{{\mathit{f}}_{\lfloor {\varphi }_{\ell }\rceil ,\ell }\}}_{\ell =1:T}\hfill \\ (\mathrm{background}\mathrm{only})\hfill & \\ {\mathcal{H}}_1:{\{{\mathit{r}}_{\lfloor {\varphi }_{\ell }\rceil ,\ell }\}}_{\ell =1:T}\hfill & =\alpha {\{{\mathit{h}}_{\lfloor {\varphi }_{\ell }\rceil }({\varphi }_{\ell })\}}_{\ell =1:T}+{\{{\mathit{f}}_{\lfloor {\varphi }_{\ell }\rceil ,\ell }\}}_{\ell =1:T},\hfill \\ (\mathrm{background}+\mathrm{source}),\hfill \end{array}\end{array}$$(2)

where α is the flux of the point source and h_⌊ϕℓ⌉(ϕ_ℓ) is a patch of the off-axis PSF extracted around the integer location ⌊ϕ_ℓ⌉, for a point source located at the subpixel location ϕ_ℓ.

Under the assumption that each of the T background patches of the collection {f_{⌊ϕℓ⌉, ℓ}}_{ℓ = 1 : T} is an independent realization, each distributed according to its local GSM model, the co-log-likelihood of the data under hypothesis ℋ₀ is:

$$\begin{array}{cc}\hfill -log\mathrm{p}({\{{\mathit{r}}_{\lfloor {\varphi }_{\ell }\rceil ,\ell }\}}_{\ell =1:T}|{\mathcal{H}}_0)=& \frac{\mathit{TK}}{2}log2\pi +{\displaystyle \sum _{\ell =1}^T}\frac{1}{2}log\det ({\widehat{\sigma }}_{\lfloor {\varphi }_{\ell }\rceil ,\ell }^2{\widehat{\mathbf{C}}}_{\lfloor {\varphi }_{\ell }\rceil })\hfill \\ \hfill & +{\displaystyle \sum _{\ell =1}^T}\frac{1}{2}{({\mathit{r}}_{\lfloor {\varphi }_{\ell }\rceil ,\ell }-{\widehat{\mathit{m}}}_{\lfloor {\varphi }_{\ell }\rceil })}^{\mathrm{t}}\hfill \\ \hfill & \times {({\widehat{\sigma }}_{\lfloor {\varphi }_{\ell }\rceil ,\ell }^2{\widehat{\mathbf{C}}}_{\lfloor {\varphi }_{\ell }\rceil })}^{-1}({\mathit{r}}_{\lfloor {\varphi }_{\ell }\rceil ,\ell }-{\widehat{\mathit{m}}}_{\lfloor {\varphi }_{\ell }\rceil }).\hfill \end{array}$$(3)

The co-log-likelihood of the data under hypothesis ℋ₁ cannot be directly evaluated because the flux α of the point source is not known beforehand. It can however be estimated in the maximum likelihood sense; see Eq. (14) in Flasseur et al. (2018a):

$$\begin{array}{c}\hfill \widehat{\alpha }={\arg \min }_{\alpha }-logp({\{{\mathit{r}}_{\lfloor {\varphi }_{\ell }\rceil ,\ell }\}}_{\ell =1:T}|{\mathcal{H}}_1,\alpha )=\frac{{\sum }_{\ell =1}^T{b}_{\ell }}{{\sum }_{\ell =1}^T{a}_{\ell }},\end{array}$$(4)

with

$$\begin{array}{c}\hfill a_{\ell }=\frac{1}{{\widehat{\sigma }}_{\lfloor {\varphi }_{\ell }\rceil ,\ell }^2}{\mathit{h}}_{\lfloor {\varphi }_{\ell }\rceil }{({\varphi }_{\ell })}^{\mathrm{t}}·{\widehat{\mathbf{C}}}_{\lfloor {\varphi }_{\ell }\rceil }^{-1}·{\mathit{h}}_{\lfloor {\varphi }_{\ell }\rceil }({\varphi }_{\ell }),\end{array}$$(5)

and

$$\begin{array}{c}\hfill b_{\ell }=\frac{1}{{\widehat{\sigma }}_{\lfloor {\varphi }_{\ell }\rceil ,\ell }^2}{\mathit{h}}_{\lfloor {\varphi }_{\ell }\rceil }{({\varphi }_{\ell })}^{\mathrm{t}}·{\widehat{\mathbf{C}}}_{\lfloor {\varphi }_{\ell }\rceil }^{-1}·({\mathit{r}}_{\lfloor {\varphi }_{\ell }\rceil ,t_{\ell }}-{\widehat{\mathit{m}}}_{\lfloor {\varphi }_{\ell }\rceil }).\end{array}$$(6)

The co-log-likelihood of the data under hypothesis ℋ₁, for a source flux $\widehat{\alpha }$, is then:

$$\begin{array}{cc}\hfill -logp({\{{\mathit{r}}_{\lfloor {\varphi }_{\ell }\rceil ,\ell }\}}_{\ell =1:T}|{\mathcal{H}}_1,\widehat{\alpha })=& \frac{\mathit{TK}}{2}log2\pi +{\displaystyle \sum _{\ell =1}^T}\frac{1}{2}log\det ({\widehat{\sigma }}_{\lfloor {\varphi }_{\ell }\rceil ,\ell }^2{\widehat{\mathbf{C}}}_{\lfloor {\varphi }_{\ell }\rceil })\hfill \\ \hfill & +{\displaystyle \sum _{\ell =1}^T}\frac{1}{2}{\mathit{u}}_{\lfloor {\varphi }_{\ell }\rceil ,\ell }^{\mathrm{t}}{({\widehat{\sigma }}_{\lfloor {\varphi }_{\ell }\rceil ,\ell }^2{\widehat{\mathbf{C}}}_{\lfloor {\varphi }_{\ell }\rceil })}^{-1}{\mathit{u}}_{\lfloor {\varphi }_{\ell }\rceil ,\ell },\hfill \end{array}$$(7)

with ${\mathit{u}}_{\lfloor {\varphi }_{\ell }\rceil ,\ell }={\mathit{r}}_{\lfloor {\varphi }_{\ell }\rceil ,\ell }-\widehat{\alpha }{\mathit{h}}_{\lfloor {\varphi }_{\ell }\rceil }({\varphi }_{\ell })-{\widehat{\mathit{m}}}_{\lfloor {\varphi }_{\ell }\rceil }$.

After some simplifications, the generalized likelihood ratio test (GLRT) is obtained; see Kay (1998) and Flasseur et al. (2018a):

$$\begin{array}{c}\hfill 2log\frac{\mathrm{p}({\{{\mathit{r}}_{\lfloor {\varphi }_{\ell }\rceil ,\ell }\}}_{\ell =1:T}|{\mathcal{H}}_1,\widehat{\alpha })}{\mathrm{p}({\{{\mathit{r}}_{\lfloor {\varphi }_{\ell }\rceil ,\ell }\}}_{\ell =1:T}|{\mathcal{H}}_0)}=\frac{({\sum }_{\ell =1}^T{b}_{\ell }{)}^2}{{\sum }_{\ell =1}^T{a}_{\ell }}\underset{{\mathcal{H}}_0}{\stackrel{{\mathcal{H}}_1}{\gtrless }}\eta .\end{array}$$(8)

Only a positive flux α makes sense for a source. The estimate $\widehat{\alpha }$ should then be replaced by the estimate ${\widehat{\alpha }}^{+}=\max (\widehat{\alpha },0)$ obtained under a positivity constraint. The GLRT can then be recast, for η ≥ 0, into the form:

$$\begin{array}{c}\hfill S/N\mathrm{test}:\frac{{\sum }_{\ell =1}^T{b}_{\ell }}{\sqrt{\sum _{\ell =1}^T{a}_{\ell }}}=\frac{\widehat{\alpha }}{{\widehat{\sigma }}_{\alpha }}\underset{{\mathcal{H}}_0}{\stackrel{{\mathcal{H}}_1}{\gtrless }}\tau ,\end{array}$$(9)

with $\tau =\sqrt{\eta }$, which corresponds to comparing the signal-to-noise ratio (S/N) of the flux estimate to a threshold; see Mugnier et al. (2009) and Flasseur et al. (2018a). The S/N test defined in (9) has the attractive property of being an affine transformation of the observations. Under our GSM model, and thanks to the normalization by the scaling factors ${\widehat{\sigma }}_{\lfloor {\varphi }_{\ell }\rceil ,\ell }$, the ratio $\widehat{\alpha }/{\widehat{\sigma }}_{\alpha }$ is distributed according to a standard normal distribution. This simplifies the setting of a detection threshold τ: τ = 5 leads to a probability of false alarms equal to 2.87 × 10⁻⁷, that is, the probability of a Gaussian random variable being greater than five standard deviations.

Figure 4 illustrates that the ratio $\widehat{\alpha }/{\widehat{\sigma }}_{\alpha }$ is indeed distributed like a standard normal distribution. Two detection maps are shown for the star HIP 72192. The left part of the figure corresponds to the results of PACO algorithm, as presented in Flasseur et al. (2018a). The corresponding detection map is stationary. When representing the empirical distribution of S/N values in the field of view (after excluding the two areas corresponding to two sources), a good match with a Gaussian distribution is obtained, albeit with a standard deviation slightly below a value of 1. The right-hand part of the figure gives the detection map produced by the extension of PACO introduced in this paper. While both maps are quite similar, it can be noted that time-dependent scaling factors lead to a spatially more stationary map and an improved match to the standard Gaussian distribution. The impact of producing maps with a variance of less than one in the absence of a source with the standard PACO algorithm is to be overly conservative: some detections may be missed at a given false alarm rate while they would be correctly detected with the robust version of PACO. We illustrate such behavior in Sect. 4.

The detection map is produced by computing the S/N for each reference location ϕ₀ on a grid. The computational complexity can be significantly reduced by pre-computing terms that can then be shared to evaluate the S/N at several locations. Such an approach, named “fast PACO” in Flasseur et al. (2018a), can also be adapted to our GSM model, provided that the background statistics are computed using the algorithm outlined in Sect. 3.1 and that the formulas are updated to include the scaling factors ${\widehat{\sigma }}_{k,\ell }$ according to Eqs. (5), (6), and (9). Compared to the original fast PACO algorithm, these changes lead to an increase of the computation time by an order of magnitude, corresponding to the typical number of iterations to reach a fixed point in the step that estimates the local background statistics (about ten iterations).

Fig. 4. Detection map computed on a SPHERE-IRDIS ADI dataset of HIP 72192. The two faint point sources in the field of view, FPS1 and FPS2, are marked by black dots. There are no other detectable sources in the field of view. Left column: standard PACO algorithm, right column: proposed extension to improve robustness. Bottom row: solid line showing the empirical distribution of the S/N values of the detection maps, excluding areas FPS1 and FPS2. The dashed line corresponds to a standard normal distribution.

3.3. Robust estimation of photometry and astrometry

In Appendix B, we prove that the GSM model considered in this paper leads to robust estimates of the background statistics and of the flux of a point source in the sense that the estimates remain bounded when a frame takes arbitrarily large values. This contrasts with other data-processing methods in ADI that generally lead to incorrect flux estimates if bad frames with large values are present in the dataset³.

Similarly to the approach in Flasseur et al. (2018a), when characterizing a source found above the detection threshold in the detection map, the source flux is to be re-estimated jointly with the background statistics in order to prevent any self-subtraction that would bias the estimation. In practice, this joint estimation is performed by alternating (i) a re-estimation of the local background statistics using the residuals ${\{{\mathit{r}}_{\lfloor {\varphi }_{\ell }\rceil ,\ell }-\widehat{\alpha }{\mathit{h}}_{\lfloor {\varphi }_{\ell }\rceil }({\varphi }_{\ell })\}}_{\ell =1:T}$, where $\widehat{\alpha }$ is the current flux estimate, and (ii) a re-estimation of the source flux with the updated background statistics by applying Eq. (4).

The accurate astrometric estimation of the source is obtained by evaluating the S/N, with the unbiased flux estimate, over a refined subpixel grid around the detected location.

Once a source has been characterized both in terms of astrometry and photometry, its contribution to the data can be subtracted and the detection map updated accordingly, similarly to source extraction with the CLEAN algorithm; see a more detailed description in Flasseur et al. (2018b).

4. Characterization of detection, and astrometric and photometric performances

The maximum achievable contrast between the host star and an off-axis point source is important information to characterize the overall performance of the instrument including the data processing. Reliable evaluations of the detection limit (expressed in contrast) typically require the injection of a fake source at different locations of the field of view and for different values of contrast. This Monte-Carlo simulation approach is computationally very costly, and therefore other methods to evaluate contrast are followed. These other methods are based on a statistical modeling and lead to predictions that are valid only to the extent that the statistical model itself is valid. It was shown in Flasseur et al. (2018b) that contrast curves produced by reference algorithms like TLOCI and KLIP were over-optimistic due to a coarse statistical model of the residues. Contrast curves obtained with PACO are closer to the performance achieved by Monte-Carlo, however we pointed out in Flasseur et al. (2018a) that they only provided a lower bound, which was not reached in practice when the statistics of the background were learnt, in the detection step, on patches that contained both the background and the source(s).

Figure 5 presents contrast curves computed on a SPHERE-IRDIS ADI dataset (at λ = 2.110 μm) of HIP 72192. Four contrast curves are reported, corresponding to TLOCI, KLIP, PACO, and our robust extension of PACO. Contrast curves of TLOCI and KLIP should be analyzed with caution because they were computed based on a threshold at 5σ that does not correspond to a probability of false alarm of 2.87 × 10⁻⁷ since the residuals are not Gaussian distributed. Due to the nonstationarity of TLOCI and KLIP S/N maps, a spatially variant threshold should be used to reach a constant false-alarm rate and compute reliable contrast values. For PACO and robust PACO, we give both the theoretical contrast curve and the results of Monte-Carlo simulations computed for five different angular separations (triangles for PACO, squares for robust PACO). The theoretical contrast curves, corresponding to fluxes α such that $\mathrm{P}(\widehat{\alpha }/{\widehat{\sigma }}_{\alpha }\ge 5)=50\%$, are obtained by computing $5{\widehat{\sigma }}_{\alpha }$, that is, $5/\sqrt{\sum _{\ell }{a}_{\ell }}$, at each point of the field of view, with intensities expressed relative to the intensity of the host star. PACO and robust PACO achieve better contrast than TLOCI and KLIP. The discrepancy between the theoretical contrast and the Monte-Carlo simulations is reduced with robust PACO compared to PACO. Moreover, the contrast is clearly improved at small separations.

Fig. 5. Contrast curves at 5σ for TLOCI, KLIP, PACO, and the proposed robust extension of PACO, for HIP 72192 (SPHERE-IRDIS, λ = 2.110 μm). The five triangles and five squares correspond to Monte-Carlo injections (see text).

Figure 6 gives the 2D contrast maps for three detection methods, on the same ADI dataset as that studied in Fig. 5. From left to right are displayed contrast maps of the standard PACO algorithm, PACO algorithm after manually removing the eight frames that display the largest fluctuations of the series, and the proposed robust extension of PACO. The gains in contrast obtained with respect to standard PACO are given on the bottom row of the figure. Manually removing the bad frames has two effects: (i) an improvement of the contrast in the area impacted by the large fluctuations of the bad frames, and (ii) a slight degradation of the contrast in the rest of the field of view due to the reduction of the size of the ADI dataset. By locally adapting the scaling factors, robust PACO improves the contrast everywhere. Even in regions that are not affected by strong temporal heterogeneity, it is beneficial to account for fluctuations of the scaling factors. As expected, the gain is the largest in the area affected by the strongest fluctuations. There, the gain in contrast reaches 80%.

Fig. 6. Comparison of contrast maps for HIP 72192 (SPHERE-IRDIS). Top row: largest contrast between the host star and a point source that would theoretically lead to a probability of detection greater than or equal to 50% at the detection threshold τ = 5. Three methods are compared: the standard PACO algorithm applied on the whole ADI stack (left), PACO applied to the subset of the ADI stack obtained by removing the eight frames with the largest fluctuations (center), and the proposed robust extension of PACO applied on the whole stack (right). The bottom row gives the gain in contrast with respect to PACO algorithm.

Standard PACO and its robust extension are also compared on the same SPHERE-IRDIS dataset as introduced in Flasseur et al. (2018a). In this dataset, 30 fake point sources were injected in addition to the two known faint point sources, at contrasts such that they could all be detected without a false alarm (see Table 3 of Flasseur et al. (2018a) for the corresponding levels of contrast). Figure 7 gives the detection maps produced by PACO and robust PACO. We refer the reader to Fig. 10 in Flasseur et al. (2018a) for a comparison with other detection methods. By including time-specific scaling factors, robust PACO improves the S/N of all sources. The largest increase occurs close to the host star where the S/N is improved by more than 100% for two point sources.

Fig. 7. Comparison of PACO (left) and the robust extension of PACO introduced in this paper (right) on the SPHERE-IRDIS ADI dataset with 30 fake point sources injected in addition to the two known faint point sources (same level of contrast as presented in Flasseur et al. 2018a). The detection maps are given together with the S/N values of the 60 first detections. The six fake sources closest to the host star are denoted by digits ① to ⑥. The two real sources are indicated in pink. The increase in S/N brought by the robust approach is given as a percentage on the detection map.

Figure 8 compares TLOCI, KLIP, standard PACO, and robust PACO algorithms on an ADI stack of HD 95086 (T = 52 temporal frames acquired under the 2015-05-05 – 095.C-0298(A) ESO program with a total apparent rotation of the field of view of 18.2°). The HD 95086 hosts an exoplanet of 5M_J in mass (HD 95086 b) discovered (Rameau et al. 2013a) and confirmed (Rameau et al. 2013b) by direct imaging with the SPHERE instrument. In addition, six known background point sources are in the field of view. Figure 8 shows that both the PACO method and its robust extension achieve better detection performance than reference algorithms; in particular, the S/N values are larger for all sources (except for the background star 2 which poses no detection problem in itself: it is so bright that it is visible directly in a single image of the ADI dataset). With PACO and robust PACO, it is possible to detect without any false alarm the seven known sources by thresholding the S/N maps at τ = 5. In comparison, S/N maps from TLOCI and KLIP present several artifacts due to the presence of outliers in the datasets, preventing an automatic analysis of these maps. The robust extension of PACO increases the detection confidence of all sources present in the field of view, especially for the ones with the smallest angular separations, such as the exoplanet HD 95086 b (S/N = 6.1 with PACO vs. 8.5 with robust PACO).

Fig. 8. Maps of S/N from TLOCI, KLIP, PACO, and the robust extension of PACO introduced in this paper on an ADI dataset of HD 95086. Seven circular insets show a zoom of the S/N maps around the known point sources. The S/N maps of these insets are evaluated on a subpixel grid (four nodes per pixel) for PACO and robust PACO. True detections are marked by straight circles, the missed detections (S/N <  5) are marked by dashed circles, and false alarms are identified by red squares.

The presence of frames of poorer quality in an ADI stack degrades the detection performance with respect to an ADI stack of constant quality. It may be useful to the astronomers who planned the observations to compare the contrast achieved given the observation conditions experienced to the contrast that would have been reached should the conditions have been as good as that of the best frames. Such a difference in contrast can be assessed by computing a “best case” contrast where ${\widehat{\sigma }}_{\lfloor {\varphi }_{\ell }\rceil ,\ell }^2$ in the computation of a_ℓ in Eq. (5) is replaced by $\underset{{\ell }^{\prime }}{min}{\widehat{\sigma }}_{\lfloor {\varphi }_{\ell }\rceil ,{\ell }^{\prime }}^2$ (i.e., the smallest scaling factor at that location ⌊ϕ_ℓ⌉). Figure 9 shows the “best case” contrast for HIP 72192. Compared to the actual contrast predicted based on the estimated time-specific scaling factors, the “best case” contrast is about 25% better in most of the field of view, and close to 70–80% in the region most impacted by large temporal fluctuations. If the aim of the observations was to characterize a faint point source that falls in that area, it might be worth considering re-observing in the hope of experiencing better conditions and getting closer to the best-case contrast.

Fig. 9. Left: contrast that would be obtained for HIP 72192 under ideal observation conditions (i.e., as good as that of the best image, during the entire ADI acquisition). Right: gain in contrast compared to the actual contrast, computed by accounting for temporal fluctuations of the scale of the residuals.

The astrometric and photometric accuracies can be characterized over the whole field of view by computing the Cramér-Rao lower bounds of the vector p of parameters (the 2D angular location ϕ₀ and the flux α) that characterize a point source. The Cramér-Rao lower bound is a good estimate of the covariance of the maximum likelihood estimator when the number of samples is large enough (Kendall et al. 1948). We follow the approach of Flasseur et al. (2018a): the Fisher information matrix I^F on the vector p is given by

$$\begin{array}{c}\hfill [{\mathbf{I}}^{\mathrm{F}}(\mathit{p}){]}_{i,j}={\displaystyle \sum _{\ell =1}^T}{\frac{1}{{\widehat{\sigma }}_{\lfloor {\varphi }_0\rceil ,\ell }^2}\frac{\partial \left[\alpha {\mathit{h}}_{\lfloor {\varphi }_0\rceil }({\varphi }_0)\right]}{\partial p_i}}^{\mathrm{t}}{\widehat{\mathbf{C}}}_{\lfloor {\varphi }_0\rceil }^{-1}\frac{\partial \left[\alpha {\mathit{h}}_{\lfloor {\varphi }_0\rceil }({\varphi }_0)\right]}{\partial p_i},\end{array}$$(10)

where i and j vary in value between 1 and 3 (p₁ and p₂ represent the two components of the 2D location ϕ₀ and p₃ corresponds to the flux α). The product α h_θk models the signal of a point source of flux α. As in Flasseur et al. (2018a), we use a continuous model of the off-axis PSF (an isotropic Gaussian can be used) to simplify the computation of the spatial derivatives. The standard deviations δ_i for each of the three parameters are obtained from the diagonal of the inverse of Fisher information matrix:

$$\begin{array}{c}\hfill [\mathit{\delta }(\mathit{p}){]}_i=\sqrt{[{\mathbf{I}}^{\mathrm{F}}{(\mathit{p})}^{-\mathbf{1}}{]}_{i,i}}.\end{array}$$(11)

Figures 10 and 11 give the astrometric Cramér-Rao lower bounds on the whole field of view obtained with PACO and robust PACO. They show that the patch weighting included in the robust approach improves the estimation accuracies on the whole field of view, especially in the presence of high stellar leakages (more than 50% improvement near the coronagraph).

Fig. 10. Theoretical accuracy (minimal standard deviation given by the Cramér-Rao lower bounds) δx0 on the estimated astrometric location x0 from PACO and its robust extension presented in this paper around the HIP 72192 star. Here, δx0 maps are multiplied by the flux α of the exoplanet (expressed in arcsec × flux). The gain (in percent) brought by the robust approach is given on the right. The six fake sources closest to the host star are denoted by digits ① to ⑥.

Fig. 11. Theoretical accuracy (minimal standard deviation given by the Cramér-Rao lower bounds) δy0 on the estimated astrometric location y0 from PACO and its robust extension presented in this paper around the HIP 72192 star. Here, δy0 maps are multiplied by the flux α of the exoplanet (expressed in arcsec × flux). The gain (in percent) brought by the robust approach is given on the right. The six fake sources closest to the host star are denoted by digits ① to ⑥.

5. Conclusion

Most data-processing techniques for point-source detection and characterization in ADI datasets are not robust: the presence of images with poor AO correction can considerably degrade their performance. It is therefore necessary to detect and eliminate these frames in a preprocessing step. In contrast, the approach described in this paper models the temporal variations of the amplitude of the background fluctuations jointly to the spatial correlations using a Gaussian scale mixture model. Since these fluctuations are accounted for, the estimations are robust to large fluctuations. Rather than discarding full images, we use a local modeling in order to spatially adapt the processing and give a very small statistical weight only to the areas subject to large fluctuations. That way, the proposed robust PACO processing makes the most of the available data.

The robust variant of the PACO algorithm introduced in this paper leads to improved detection performances, in particular at close separations where the stellar residuals dominate. Interestingly, it is also possible to estimate the achievable contrast not only by taking into account the actual image quality of the ADI dataset, but also the contrast that would have been reached if the AO correction had been as good as the best correction in the ADI time series. This information can be highly valuable in order to plan subsequent re-observations.

Coronagraphic observations obtained with integral field spectrographs could also benefit from a modeling of the spatial and temporal fluctuations. We will investigate in a future work how to include a modeling of the spectral correlations for these instruments.

Acknowledgments

We thank A. Boccaletti (LESIA, Paris, France) who provided the transmission of the SPHERE coronagraphs. We also thank the anonymous Referee for his careful reading of the manuscript as well as his insightful comments and suggestions. This work has made use of the SPHERE Data Centre, jointly operated by OSUG/IPAG (Grenoble, France), PYTHEAS/LAM/CESAM (Marseille, France), OCA/Lagrange (Nice, France), Observatoire de Paris/LESIA (Paris, France), and Observatoire de Lyon/CRAL (Lyon, France).

References

Appendix A: Estimation of the local statistics of the background

Let {r_θk, tℓ}_{ℓ = 1 : T} be a collection of 2D patches centered at location θ_k. Under our model, the ℓ-th patch of the collection is distributed under ℋ₀ according to a normal distribution with mean m_k and spatial covariance ${\mathbf{C}}_{k,\mathcal{l}}={\sigma }_{k,\mathcal{l}}^2{\mathbf{C}}_k$. In this Appendix, we derive the expression of the maximum likelihood estimators ${\widehat{\mathit{m}}}_k$, ${\widehat{\sigma }}_{k,\ell }^2$ and ${\widehat{\mathbf{S}}}_k$ of m_k, ${\sigma }_{k,\mathcal{l}}^2$ and C_k, respectively. In the paper, the notation ${\widehat{\mathbf{C}}}_k$ is kept to denote the shrinkage estimator obtained by reducing the value of off-diagonal elements of the maximum likelihood estimator ${\widehat{\mathbf{S}}}_k$.

Under the assumption that each patch is an independent realization of our model of the background, the co-log-likelihood ℒ of the collection is given by:

$$\begin{array}{cc}\hfill \mathcal{L}=& -log\mathrm{p}({\{{\mathit{r}}_{{\theta }_k,t_{\ell }}\}}_{\ell =1:T}|{\mathit{m}}_k,{\{{\sigma }_{k,\ell }^2\}}_{\ell =1:T},{\mathbf{C}}_k)\hfill \\ \hfill & \propto {\displaystyle \sum _{\ell =1}^T}log\det \left({\sigma }_{k,\ell }^2{\mathbf{C}}_k\right)+{\displaystyle \sum _{\ell =1}^T}\frac{1}{{\sigma }_{k,\ell }^2}{\overline{\mathit{r}}}_{{\theta }_k,t_{\ell }}^{\mathrm{t}}{\mathbf{C}}_k^{-1}{\overline{\mathit{r}}}_{{\theta }_k,t_{\ell }}+\mathrm{const}.,\hfill \end{array}$$(A.1)

with ${\overline{\mathit{r}}}_{{\theta }_k,t_{\ell }}={\mathit{r}}_{{\theta }_k,t_{\ell }}-{\mathit{m}}_k$. The maximum likelihood estimators correspond to the location of the minimum of (A.1), where the gradient of the co-log-likelihood is equal to zero.

From the condition ${\frac{\partial \mathcal{L}}{\partial {\sigma }_{k,\ell }^2}|}_{{\sigma }_{k,\ell }^2={\widehat{\sigma }}_{k,\ell }^2}=0$, we get:

$$\begin{array}{c}\hfill \frac{K}{{\widehat{\sigma }}_{k,\ell }^2}-\frac{1}{{\widehat{\sigma }}_{k,\ell }^4}{\overline{\mathit{r}}}_{{\theta }_k,t_{\ell }}^{\mathrm{t}}{\mathbf{C}}_k^{-1}{\overline{\mathit{r}}}_{{\theta }_k,t_{\ell }}=0,\end{array}$$

which leads to

$$\begin{array}{c}\hfill {\widehat{\sigma }}_{k,\ell }^2=\frac{1}{K}{\overline{\mathit{r}}}_{{\theta }_k,t_{\ell }}^{\mathrm{t}}{\mathbf{C}}_k^{-1}{\overline{\mathit{r}}}_{{\theta }_k,t_{\ell }}.\end{array}$$(A.2)

The variance at time t_ℓ around pixel k is therefore estimated by computing the sample variance of a spatially whitened version of the ℓ-th patch.

From the condition ${{\mathbf{\nabla }}_{{\mathit{m}}_k}\mathcal{L}|}_{{\mathit{m}}_k={\widehat{\mathit{m}}}_k}=\mathbf{0}$, we get:

$$\begin{array}{c}\hfill -2{\displaystyle \sum _{\ell =1}^T}\frac{1}{{\sigma }_{k,\ell }^2}{\mathbf{C}}_k^{-1}({\mathit{r}}_{{\theta }_k,t_{\ell }}-{\widehat{\mathit{m}}}_k)=\mathbf{0}.\end{array}$$

Since ${\mathbf{C}}_k^{-1}$ is not singular, we obtain:

$$\begin{array}{c}\hfill {\widehat{\mathit{m}}}_k=\frac{{\sum }_{\ell =1}^T\frac{1}{{\sigma }_{k,\ell }^2}{\mathit{r}}_{{\theta }_k,t_{\ell }}}{{\sum }_{\ell =1}^T\frac{1}{{\sigma }_{k,\ell }^2}}·\end{array}$$(A.3)

The mean background in patches around the k-th pixel corresponds to a weighted average where a patch at time t_ℓ with large variance is given a small weight.

From the condition ${{\mathbf{\nabla }}_{{\mathbf{C}}_k}\mathcal{L}|}_{{\mathbf{C}}_k={\widehat{\mathbf{S}}}_k}=\mathbf{0}$, we derive:

$$\begin{array}{c}\hfill T·{\widehat{\mathbf{S}}}_k^{-1}-{\widehat{\mathbf{S}}}_k^{-1}\left({\displaystyle \sum _{\ell =1}^T}\frac{1}{{\sigma }_{k,\ell }^2}{\overline{\mathit{r}}}_{{\theta }_k,t_{\ell }}{\overline{\mathit{r}}}_{{\theta }_k,t_{\ell }}^{\mathrm{t}}\right){\widehat{\mathbf{S}}}_k^{-1}=\mathbf{0},\end{array}$$

which gives:

$$\begin{array}{c}\hfill {\widehat{\mathbf{S}}}_k=\frac{1}{T}{\displaystyle \sum _{\ell =1}^T}\frac{1}{{\sigma }_{k,\ell }^2}{\overline{\mathit{r}}}_{{\theta }_k,t_{\ell }}{\overline{\mathit{r}}}_{{\theta }_k,t_{\ell }}^{\mathrm{t}}.\end{array}$$(A.4)

The spatial covariance estimator corresponds to the sample covariance upon a proper rescaling by a factor 1/σ_k, ℓ of residual patches ${\overline{\mathit{r}}}_{{\theta }_k,t_{\ell }}$.

Appendix B: Proof of the robustness of the estimators

In this appendix, we show that our estimators are robust in the sense that, should the patch r_{θk, tℓ^⋆} be replaced by a scaled version ${\mathit{r}}_{{\theta }_k,t_{{\mathcal{l}}^{\star }}}^{\prime }=\eta \cdot {\mathit{r}}_{{\theta }_k,t_{{\mathcal{l}}^{\star }}}$ with an arbitrarily large factor η, at a given frame ℓ^⋆, the estimators of the background statistics ${\widehat{\mathit{m}}}_k^{\prime }$ and ${\widehat{\mathbf{S}}}_k^{\prime }$ as well as the flux estimator ${\widehat{\alpha }}^{\prime }$ computed on the collection with the scaled patch remain bounded. In fact, when η grows to infinity, we show that the patch ${\mathit{r}}_{{\theta }_k,t_{{\mathcal{l}}^{\star }}}^{\prime }$ gets completely discarded in the estimation of the flux.

Since the estimators ${\widehat{\mathit{m}}}_k$, ${\widehat{\mathbf{C}}}_k$ and ${\widehat{\sigma }}_{k,\ell }^2$ are defined implicitly through a fixed-point procedure, directly proving the robustness is difficult. Instead, we show that when η approaches infinity, the following estimates correspond to a fixed point:

$$\begin{array}{c}\hfill \{\begin{array}{cc}{\widehat{\sigma }}_{k,\ell }^{2\prime }={\stackrel{\sim }{\sigma }}_{k,\ell }^2\hfill & \mathrm{when}\ell \ne {\ell }^{\star }\hfill \\ {\widehat{\sigma }}_{k,{\ell }^{\star }}^{2\prime }\to \infty \hfill \\ {\widehat{\mathit{m}}}_k^{\prime }={\stackrel{\sim }{\mathit{m}}}_k\hfill \\ {\widehat{\mathbf{S}}}_k^{\prime }=\frac{T-1}{T}{\stackrel{\sim }{\mathbf{S}}}_k+\frac{K}{T}\frac{{\mathit{r}}_{{\theta }_k,t_{{\ell }^{\star }}}{\mathit{r}}_{{\theta }_k,t_{{\ell }^{\star }}}^{\mathrm{t}}}{{\mathit{r}}_{{\theta }_k,t_{{\ell }^{\star }}}^{\mathrm{t}}{({{\widehat{\mathbf{S}}}_k}^{\prime })}^{-1}{\mathit{r}}_{{\theta }_k,t_{{\ell }^{\star }}}},\hfill \end{array}\end{array}$$(B.1)

where ${\stackrel{\sim }{\sigma }}_{k,\ell }^2$, ${\stackrel{\sim }{\mathit{m}}}_k$ and ${\stackrel{\sim }{\mathbf{S}}}_k$ correspond to the estimators computed on the collection {r_θk, tℓ}_{ℓ ≠ ℓ^⋆} of the patches with the ℓ^⋆-th frame removed.

Starting with the initial values given in (B.1), the application of Eq. (A.2) leaves the estimated variances ${\widehat{\sigma }}_{k,\ell }^{2\prime }$ unchanged for all ℓ different from ℓ^⋆. For ℓ^⋆, ${\widehat{\sigma }}_{k,\ell }^{2\prime }\underset{\eta \to \infty }{\sim }{\eta }^2{\widehat{\sigma }}_{k,\ell }^2$ which tends to infinity when η approaches infinity.

From Eq. (A.3), we see that the patch ${\mathit{r}}_{{\theta }_k,t_{{\mathcal{l}}^{\star }}}^{\prime }$ with associated infinite variance ${\widehat{\sigma }}_{k,\ell }^{2\prime }$ has zero weight: $\frac{1}{{\widehat{\sigma }}_{k,{\ell }^{\star }}^{2\prime }}{\mathit{r}}_{{\theta }_k,t_{{\ell }^{\star }}}^{\prime }\underset{\eta \to \infty }{\sim }\frac{1}{\eta {\widehat{\sigma }}_{k,\ell }^2}{\mathit{r}}_{{\theta }_k,t_{{\ell }^{\star }}}$, so that $\underset{\eta \to \infty }{lim}\frac{1}{{\widehat{\sigma }}_{k,{\ell }^{\star }}^{2\prime }}{\mathit{r}}_{{\theta }_k,t_{{\ell }^{\star }}}^{\prime }=0$. The mean patch then corresponds to that obtained on the collection of patches with the ℓ^⋆-th patch excluded.

The approximations ${\widehat{\sigma }}_{k,{\ell }^{\star }}^{2\prime }\underset{\eta \to \infty }{\sim }\frac{{\eta }^2}{K}{\mathit{r}}_{{\theta }_k,t_{{\ell }^{\star }}}^{\mathrm{t}}{({\widehat{\mathbf{S}}}_k^{\prime })}^{-1}{\mathit{r}}_{{\theta }_k,t_{{\ell }^{\star }}}$ and ${\overline{\mathit{r}}}_{{\theta }_k,{\ell }^{\star }}^{\prime }\underset{\eta \to \infty }{\sim }\eta {\mathit{r}}_{{\theta }_k,{\ell }^{\star }}$ together with the update of the sample covariance estimate (A.4) lead to ${\widehat{\mathbf{S}}}_k^{\prime }=\frac{T-1}{T}{\stackrel{\sim }{\mathbf{S}}}_k+\frac{K}{T}\frac{{\mathit{r}}_{{\theta }_k,t_{{\ell }^{\star }}}{\mathit{r}}_{{\theta }_k,t_{{\ell }^{\star }}}^{\mathrm{t}}}{{\mathit{r}}_{{\theta }_k,t_{{\ell }^{\star }}}^{\mathrm{t}}{({\widehat{\mathbf{S}}}_k^{\prime })}^{-1}{\mathit{r}}_{{\theta }_k,t_{{\ell }^{\star }}}}$.

In conclusion, the set of parameters defined in (B.1) is a fixed-point. Moreover, neither the mean background patch ${\widehat{m}}_k^{\prime }$ nor the spatial covariance ${\widehat{S}}_k^{\prime }$ depend on η (i.e., they both remain bounded when η → ∞). It remains to show the robustness of the flux estimator $\widehat{\alpha }$ of a source.

If the scaled patch is not in the collection of all patches that contain the PSF of the source, then its influence is limited since we have just shown that the background statistics are robustly estimated. In the case of the superimposition of a scaled background and of the source of interest (as is the case when a poor correction of the AO leads to large stellar leakages at some locations of the field of view), the estimated flux becomes:

$$\begin{array}{c}\hfill \widehat{\alpha }({\varphi }_{t_{\ell }})=\frac{{\displaystyle \sum _{\ell =1}^T}\frac{1}{{\widehat{\sigma }}_{k,\ell }^{2\prime }}{\mathit{h}}_{\lfloor {\varphi }_{t_{\ell }}\rceil }{({\varphi }_{t_{\ell }})}^{\mathrm{t}}·{\widehat{\mathbf{C}}}_{\lfloor {\varphi }_{t_{\ell }}\rceil }^{-1}·({\mathit{r}}_{\lfloor {\varphi }_{t_{\ell }}\rceil ,t_{\ell }}^{\mathbf{\prime }}-{\widehat{\mathit{m}}}_{\lfloor {\varphi }_{t_{\ell }}\rceil })}{{\displaystyle \sum _{\ell =1}^T}\frac{1}{{\widehat{\sigma }}_{k,\ell }^{2\prime }}{\mathit{h}}_{\lfloor {\varphi }_{t_{\ell }}\rceil }{({\varphi }_{t_{\ell }})}^{\mathrm{t}}·{\widehat{\mathbf{C}}}_{\lfloor {\varphi }_{t_{\ell }}\rceil }^{-1}·{\mathit{h}}_{\lfloor {\varphi }_{t_{\ell }}\rceil }({\varphi }_{t_{\ell }})}·\end{array}$$(B.2)

Therefore, when η becomes large, the term $({\mathit{r}}_{\lfloor {\varphi }_{t_{{\ell }^{\star }}}\rceil ,t_{{\ell }^{\star }}}^{\prime }-{\widehat{\mathit{m}}}_{\lfloor {\varphi }_{t_{{\ell }^{\star }}}\rceil })/{\widehat{\sigma }}_{k,{\ell }^{\star }}^{2\prime }$ tends to zero and we obtain:

$$\begin{array}{c}\hfill \widehat{\alpha }({\varphi }_{t_{\ell }})\approx \frac{\sum _{\ell \ne {\ell }^{\star }}\frac{1}{{\widehat{\sigma }}_{k,\ell }^{2\prime }}{\mathit{h}}_{\lfloor {\varphi }_{t_{\ell }}\rceil }{({\varphi }_{t_{\ell }})}^{\mathrm{t}}·{\widehat{\mathbf{C}}}_{\lfloor {\varphi }_{t_{\ell }}\rceil }^{-1}·({\mathit{r}}_{\lfloor {\varphi }_{t_{\ell }}\rceil ,t_{\ell }}^{\mathbf{\prime }}-{\widehat{\mathit{m}}}_{\lfloor {\varphi }_{t_{\ell }}\rceil })}{\sum _{\ell \ne {\ell }^{\star }}\frac{1}{{\widehat{\sigma }}_{k,\ell }^{2\prime }}{\mathit{h}}_{\lfloor {\varphi }_{t_{\ell }}\rceil }{({\varphi }_{t_{\ell }})}^{\mathrm{t}}·{\widehat{\mathbf{C}}}_{\lfloor {\varphi }_{t_{\ell }}\rceil }^{-1}·{\mathit{h}}_{\lfloor {\varphi }_{t_{\ell }}\rceil }({\varphi }_{t_{\ell }})},\end{array}$$(B.3)

which indicates that the scaled frame ${\mathit{r}}_{\lfloor {\varphi }_{t_{{\mathcal{l}}^{\star }}}\rceil ,t_{{\mathcal{l}}^{\star }}}^{\prime }$ tends to be completely discarded when the scaling factor η is large. Estimation of the photometry and of the astrometry is therefore robust to the presence of a frame with an arbitrarily large background.

All Figures

	Fig. 1. Central region of SPHERE-IRDIS images of HIP 72192: panel a: measured intensity for five selected frames showing strong temporal fluctuations; panel b: map of the estimated temporal scaling parameters ${\widehat{\sigma }}_{k,\ell }$ for the matching frames: ℓ = 1, 3, 8, 13 and 14.
In the text

Fig. 2. Temporal evolution and empirical distribution of the intensities within a patch: (a)–(b) in patches extracted at a location θk close to the coronagraph; (c)–(d) in patches extracted at a location θk farther from the coronagraph; (e) in all the patches from the field of view. Rows a and c: values within the patches for ten specific frames: t1, t3, t8, t13, t14, t25, t40, t54, t63 and t69. Rows b, d, and e: empirical histograms computed over all frames t1 to t69 and, in dashed line, a standard Gaussian. The first column corresponds to centered patches ${\{{\mathit{r}}_{{\theta }_k,t_{\ell }}-{\widehat{\mathit{m}}}_k\}}_{\ell =1:T}$. The second column corresponds to centered patches, after whitening the spatial correlations. The last column corresponds to centered patches that have been both whitened for the spatial correlations and equalized with the temporal scaling factors.

In the text

Fig. 3. Illustration of the convergence of Algorithm 1: 1000 random initializations of the weights ${\{{\widehat{w}}_{\ell }\}}_{\ell =1:T}$ converge to the same solution as that obtained by initializing with constant weights. Convergence at three locations θk depicted on the field of view is shown, covering three different cases: small dispersion of the value of scaling parameters (left), medium dispersion (center), and large dispersion close to the coronagraph (right).

In the text

Fig. 4. Detection map computed on a SPHERE-IRDIS ADI dataset of HIP 72192. The two faint point sources in the field of view, FPS1 and FPS2, are marked by black dots. There are no other detectable sources in the field of view. Left column: standard PACO algorithm, right column: proposed extension to improve robustness. Bottom row: solid line showing the empirical distribution of the S/N values of the detection maps, excluding areas FPS1 and FPS2. The dashed line corresponds to a standard normal distribution.

In the text

	Fig. 5. Contrast curves at 5σ for TLOCI, KLIP, PACO, and the proposed robust extension of PACO, for HIP 72192 (SPHERE-IRDIS, λ = 2.110 μm). The five triangles and five squares correspond to Monte-Carlo injections (see text).
In the text

Fig. 6. Comparison of contrast maps for HIP 72192 (SPHERE-IRDIS). Top row: largest contrast between the host star and a point source that would theoretically lead to a probability of detection greater than or equal to 50% at the detection threshold τ = 5. Three methods are compared: the standard PACO algorithm applied on the whole ADI stack (left), PACO applied to the subset of the ADI stack obtained by removing the eight frames with the largest fluctuations (center), and the proposed robust extension of PACO applied on the whole stack (right). The bottom row gives the gain in contrast with respect to PACO algorithm.

In the text

Fig. 7. Comparison of PACO (left) and the robust extension of PACO introduced in this paper (right) on the SPHERE-IRDIS ADI dataset with 30 fake point sources injected in addition to the two known faint point sources (same level of contrast as presented in Flasseur et al. 2018a). The detection maps are given together with the S/N values of the 60 first detections. The six fake sources closest to the host star are denoted by digits ① to ⑥. The two real sources are indicated in pink. The increase in S/N brought by the robust approach is given as a percentage on the detection map.

In the text

Fig. 8. Maps of S/N from TLOCI, KLIP, PACO, and the robust extension of PACO introduced in this paper on an ADI dataset of HD 95086. Seven circular insets show a zoom of the S/N maps around the known point sources. The S/N maps of these insets are evaluated on a subpixel grid (four nodes per pixel) for PACO and robust PACO. True detections are marked by straight circles, the missed detections (S/N <  5) are marked by dashed circles, and false alarms are identified by red squares.

In the text

	Fig. 9. Left: contrast that would be obtained for HIP 72192 under ideal observation conditions (i.e., as good as that of the best image, during the entire ADI acquisition). Right: gain in contrast compared to the actual contrast, computed by accounting for temporal fluctuations of the scale of the residuals.
In the text

Fig. 10. Theoretical accuracy (minimal standard deviation given by the Cramér-Rao lower bounds) δx0 on the estimated astrometric location x0 from PACO and its robust extension presented in this paper around the HIP 72192 star. Here, δx0 maps are multiplied by the flux α of the exoplanet (expressed in arcsec × flux). The gain (in percent) brought by the robust approach is given on the right. The six fake sources closest to the host star are denoted by digits ① to ⑥.

In the text

Fig. 11. Theoretical accuracy (minimal standard deviation given by the Cramér-Rao lower bounds) δy0 on the estimated astrometric location y0 from PACO and its robust extension presented in this paper around the HIP 72192 star. Here, δy0 maps are multiplied by the flux α of the exoplanet (expressed in arcsec × flux). The gain (in percent) brought by the robust approach is given on the right. The six fake sources closest to the host star are denoted by digits ① to ⑥.

In the text