Robustness to bad frames in angular differential imaging: a local weighting approach

Context. The detection of exoplanets by direct imaging is very challenging. It requires an extreme adaptive-optics (AO) system and a coronagraph as well as suitable observing strategies. In angular differential imaging, the signal-to-noise ratio is improved by combining several observations. Aims. Due to the evolution of the observation conditions and of the AO correction, the quality of the observations may vary significantly during the observing sequence. It is common practice to reject images of comparatively poor quality. We aim to decipher when this selection should be performed and what its impact on detection performance is. Methods. Rather than discarding a full image, we study the local fluctuations of the signal at each frame and derive weighting maps for each frame. These fluctuations are modeled locally directly from the data through the spatio-temporal covariance of small image patches. The weights derived from the temporal variances can be used to improve the robustness of the detection step and reduce estimation errors of both the astrometry and photometry. The impact of bad frames can be analyzed by statistically characterizing the detection and estimation performance. Results. When used together with a modeling of the spatial covariances (PACO algorithm), these weights improve the robustness of the detection method. Conclusions. The spatio-temporal modeling of the background fluctuations provides a way to exploit all acquired frames. In the case of bad frames, areas with larger fluctuations are discarded by a weighting strategy and do not corrupt the detection map or the astrometric and photometric estimations. Other areas of better quality are preserved and are included to detect and characterize sources.


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. 2003Allard et al. , 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. 2013Marois et al. , 2014Soummer et al. 2012;Smith et al. 2009;Cantalloube 2016;Gonzalez et al. 2016Gonzalez et al. , 2018Flasseur 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 5 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 spatiotemporal 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(Marois et al. , 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.

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 2 .
The fluctuations are spatially structured and their magnitude varies from one image to the other in the ADI sequence.
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 N(m k , C k ), with spatially variant mean m k ∈ R K and covariance C k ∈ R 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): where σ k, > 0 is a scalar (random or deterministic) variable that acts as a scaling parameter and the vectors u k, ∼ N(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 N(m k , σ 2 k, C k ). This corresponds to considering parameter σ k, as a deterministic nuisance parameter; see Pascal et al. (2008). 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 {r θ k ,t − m k } =1:T is displayed in the left column. Based on a local estimate C k 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: t 1 , t 3 , t 8 , t 13 , t 14 , t 25 , t 40 , t 54 , t 63 and t 69 . Rows b, d, and e: empirical histograms computed over all frames t 1 to t 69 and, in dashed line, a standard Gaussian. The first column corresponds to centered patches {r θ k ,t − m k } =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. of the spatial covariance matrix, the distribution of spatially whitened and centered patches is then displayed on the central column. The Cholesky factor 1 L k (i.e., the lower-triangular matrix such that L k L t k = C −1 k ) is used to obtain the collection of spatially whitened patches { L t k (r θ k ,t − m k )} =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 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 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 timespecific 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).

Estimation of the statistics of the background
Estimating, at each spatial location, the mean background m, the spatial covariance C, and the scale parameters σ 1 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 T = 1/ σ 2 2 / 1/ σ 4 that corresponds to the variance reduction reached when performing weighted means with the weights 1/ σ 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 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 Algorithm 1: Local background statistics estimation Input: {r 1 , . . . , r T } (stack of T spatial patches) (each patch has K pixels) Output: m (mean patch) Step 1: Estimate scale parameters Step 2: Update temporal weights Step 3: Update the mean patch m ← T =1 w · r (weighted mean) Step 4: Update the covariance Algorithm 2: Estimate scale parameters Input: {r 1 , . . . , r T } (stack of T spatial patches) Input: m (estimated mean patch) (close to the coronagraph). The convergence of Algorithm 1 is assessed based on 1000 random initializations of the weights w (following a uniform distribution on [0,1]). The vector of weights at convergence when starting from a constant vector of weights (∀ , w = 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.

Robust computation of a detection map
The Algorithm 3: Shrinkage covariance estimator 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 φ 0 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 φ 0 . 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 m φ , spatial covariance C φ and scaling factor 2 σ φ , have been estimated for all locations φ corresponding to the trajectory of a hypothetical point source with reference location φ 0 , the likelihood of two hypotheses can be compared: 2 Only the scaling factor at time t is necessary for the detection.
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-loglikelihood of the data under hypothesis H 0 is: The co-log-likelihood of the data under hypothesis H 1 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): The co-log-likelihood of the data under hypothesis H 1 , for a source flux α, is then: with u φ , = r φ , − αh φ (φ ) − m φ . After some simplifications, the generalized likelihood ratio test (GLRT) is obtained; see Kay (1998) and Flasseur et al. (2018a): Only a positive flux α makes sense for a source. The estimate α should then be replaced by the estimate α + = max( α, 0) obtained under a positivity constraint. The GLRT can then be recast, for η ≥ 0, into the form: with τ = √ η, which corresponds to comparing the signalto-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 σ φ , , the ratio α/ σ α 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 −7 , that is, the probability of a Gaussian random variable being greater than five standard deviations. Figure 4 illustrates that the ratio α/ σ α 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 φ 0 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 σ k, 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).

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 3 .
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 selfsubtraction 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 {r φ , − αh φ (φ )} =1:T , where α is the current flux estimate, and (ii) a re-estimation of the source flux with the updated background statistics by applying Eq. (4). 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.
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).

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 A2, page 7 of 13 A&A 634, A2 (2020) 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 −7 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 P( α/ σ α ≥ 5) = 50%, are obtained by computing 5 σ α , that is, 5/ a , 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. 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%.
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. 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   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. 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).
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 σ 2 φ , in the computation of a in Eq. (5) is replaced by min σ 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 A2, page 9 of 13 A&A 634, A2 (2020) 50 25 0 Fig. 10. Theoretical accuracy (minimal standard deviation given by the Cramér-Rao lower bounds) δ x 0 on the estimated astrometric location x 0 from PACO and its robust extension presented in this paper around the HIP 72192 star. Here, δ x 0 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 1 to 6 .  Fig. 11. Theoretical accuracy (minimal standard deviation given by the Cramér-Rao lower bounds) δ y 0 on the estimated astrometric location y 0 from PACO and its robust extension presented in this paper around the HIP 72192 star. Here, δ y 0 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 1 to 6 . area, it might be worth considering re-observing in the hope of experiencing better conditions and getting closer to the best-case contrast.
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 φ 0 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 where i and j vary in value between 1 and 3 (p 1 and p 2 represent the two components of the 2D location φ 0 and p 3 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: 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).

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.
Therefore, when η becomes large, the term r φ t ,t − m φ t /σ 2 k, tends to zero and we obtain: which indicates that the scaled frame r φ t ,t 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.