© ESO 2021

1. Introduction

High-contrast imaging (HCI) is a method aimed at disentangling the faint signal of an exoplanet from the much brighter signal of its host star. It relies on large telescopes, adaptive optics, coronagraphs, and sophisticated data processing. While state-of-the-art HCI instruments provide near-infrared images with a high level of sensitivity and angular resolution, they are not immune to residual aberrations. These aberrations originate from time-averaged uncorrected atmospheric turbulence or from the optical train of the telescope and instrument. The resulting high and low spatial frequency structures can be problematic when searching for exoplanetary candidates since some of them, which are known as quasi-static speckles, share the same properties. Quasi-static speckles are often correlated across the exposures of a typical observation (Bloemhof et al. 2001), which leads to their variation either being not slow enough to be subtracted through instruments calibration nor fast enough to be time-averaged.

Data processing represents a key ingredient for reducing the intensity of these speckles and further increasing the sensitivity of HCI. Different post-processing techniques, coupled with observing strategies, have been proposed in recent decades. Spectral differential imaging (SDI, Marois et al. 2000; Sparks & Ford 2002) and angular differential imaging (ADI, Marois et al. 2006) are the most commonly used observing strategies in dealing with quasi-static speckles. Most ADI-based methods rely on the subtraction of a reference point spread function (PSF), which models the speckle field based on the original set of images. The resulting residual images are then co-aligned and combined, cancelling part of the residual noise while increasing the intensity of the potential exoplanet signal.

The most common subtraction methods include: ADI median-subtraction, which subtracts the median of a set of reference frames from the original set of images, the locally optimised combination of images (LOCI, Lafreniere et al. 2007), which uses a least-squares approach to optimally estimate the speckle field, and the principal component analysis (PCA/KLIP, Soummer et al. 2012; Amara & Quanz 2012), which computes the high-order singular modes of the speckle field. More recently, other speckle subtraction methods, such as non-negative matrix factorisation (NMF, Ren et al. 2018) and the local low-rank plus sparse plus Gaussian decomposition (LLSG, Gomez Gonzalez et al. 2016) have been proposed. These methods generally rely on signal-to-noise (S/N) maps to detect planetary candidates via the definition of a S/N threshold (e.g., Mawet et al. 2014; Bottom et al. 2017; Pairet et al. 2019). Other methods such as ANDROMEDA (Cantalloube et al. 2015) or KLIP-FMMF (Pueyo 2016; Ruffio et al. 2017) use a forward model of the point source to identify the planetary signal in the residual images via the maximum likelihood. These approaches aim to take into account the distortion of the point source caused by the speckle field subtraction.

A new approach, known as a regime-switching model (RSM) map (Dahlqvist et al. 2020), was recently proposed to better take advantage of the numerous PSF subtraction techniques that have been developed in the past decade. The RSM map can be seen as a multi-ADI alternative to the estimation of S/N map as it allows for the creation of a probability map based on several cubes of residuals obtained with different ADI-based PSF subtraction techniques. In contrast with most S/N map based techniques, the RSM map does not make the assumption of a Gaussian and white distribution for the speckle residuals after subtraction. It was indeed demonstrated that the speckle residual distribution is often closer to a Laplacian distribution, especially at small separations (Pairet et al. 2019; Dahlqvist et al. 2020). The RSM map takes into account the radial evolution of the speckle noise distribution and its dependence on the instrument, switching between Gaussian and Laplacian distributions. Dahlqvist et al. (2020) successfully tested this approach on a subset of PSF subtraction techniques (annular PCA, NMF, and LLSG).

The goal of this paper is to further develop the RSM approach by considering a larger set of PSF subtraction techniques that includes the forward modelling of the point source. Indeed, as for ANDROMEDA and KLIP-FMMF, the RSM map relies on a matched filter to infer the existence of planetary candidates in residuals images. However, the initial version of the algorithm, summarised in Sect. 2, uses solely an off-axis PSF for the detection. Forward modelling could significantly improve the sensitivity of the algorithm to faint companions by taking into account the distortion generated by the speckle field subtraction. We propose a method that would rely on the KLIP forward model (KLIP-FM) developed by Pueyo (2016) as well as a forward model version of LOCI for investigating the added value of forward-modelled point sources. Section 3 is devoted to the development of the two forward-model versions of the RSM algorithm, while Sect. 4 provides a performance assessment of these versions.

As mentioned above, the RSM map can accommodate several PSF subtraction techniques to generate a final probability map. That raises questions regarding the selection of the optimal set of techniques to reach the highest sensitivity as well as its dependence on the HCI instrument and on the radial distance. We compare, in Sect. 5, the performance of several set of techniques via receiver operating characteristic (ROC) curves and investigate the impact of the considered instruments on this selection, considering three state-of-the-art HCI instruments: NACO, SPHERE, and LMIRCam.

Finally, in Sect. 6, we propose an improved method for the probability estimation, relying on a forward-backward approach that allows the use of both past and future observations within the cube of residuals to generate the RSM map. The original RSM map uses a simple forward approach, which considers only past observations to build up the probabilities. We compare the performance of both approaches with standard S/N maps through the use of ROC curves and contrast curves. In Sect. 7, we present a new framework developed to compute these contrast curves, as it is not possible when dealing with probability maps to rely on the standard procedure used for S/N maps.

2. The RSM map

The RSM approach (Dahlqvist et al. 2020) relies on a two-state Markov chain to model the pixel intensity evolution inside the de-rotated cube of residuals generated by an ADI-based PSF subtraction technique. The RSM is applied annulus-wise to account for the radial evolution of the residual noise properties. A residuals time series, x_ia, is built for each annulus, a, by vectorising the selected set of patches along the time axis and then the spatial axis. The index i_a ∈ {1, …, T × L_a} provides the position of the considered patches within the cube of residuals, where T and L_a are respectively the number of frames in the cube and the number of pixels in a given annulus of radius, a. The time series is then described by a set of two equations accounting for the two considered regimes: the first regime, S_ia = 0, where x_ia is described by a residual noise following the statistics of the quasi-static speckle residuals; and a second regime, S_ia = 1, where x_ia is described by both the residual noise and the planetary signal model¹. Thus,

$$\begin{array}{c}\hfill {\mathit{x}}_{i_a}=\mu +\beta F_{i_a}\mathit{m}+{\mathit{\epsilon }}_{\mathit{s}\mathbf{,}{\mathit{i}}_{\mathit{a}}}=\{\begin{array}{cc}\mu +{\mathit{\epsilon }}_{\mathbf{0}\mathbf{,}{\mathit{i}}_{\mathit{a}}}\hfill & \text{if}{S}_{i_a}=0\hfill \\ \mu +\beta \mathit{m}+{\mathit{\epsilon }}_{\mathbf{1}\mathbf{,}{\mathit{i}}_{\mathit{a}}}\hfill & \text{if}{S}_{i_a}=1\hfill \end{array},\end{array}$$(1)

where β and m provide, respectively, the flux and a model of the planetary signal, μ is the mean of the quasi-static speckle residuals, and ε_s, ia is their time and space varying part that is characterised by the quasi-static speckle residuals statistics (see Appendix A for a summary of all the variables used for the RSM map computation). The parameter F_ia is a realisation of a two-state Markov chain that provides a short-term memory function to the model. The parameter F_ia being the realisation of a stochastic process, the time series x_ia is described via a probability-weighted sum of the values generated by Eq. (1).

The key feature of the RSM is the fact that it provides for every pixel in the selected annulus, a, a probability for their being in each regime. The RSM map is constructed based on the probabilities of being in the planetary regime, that is, S_ia = 1. Relying on a two-state Markov chain, the probability, ξ_s, ia, of being in regime s at index i_a will depend on the probability at the previous step, i_a − 1, on the likelihood of being currently in a given regime, η_s, ia, and on the transition probability between the regimes given by the matrix, p_q, s. The resulting probability of being in the regime, s, at index, i_a, is given by:

$$\begin{array}{c}\hfill {\xi }_{s,i_a}={\displaystyle \sum _{q=0}^1\frac{{\eta }_{s,i_a}{p}_{q,s}{\xi }_{1,i_a-1}}{{\sum }_{q=0}^1{\sum }_{s=0}^1{\eta }_{s,i_a}{p}_{q,s}{\xi }_{q,i_a-1}},}\end{array}$$(2)

where $\sum _{q=0}^1\sum _{s=0}^1{\eta }_{s,i_a}{p}_{q,s}{\xi }_{q,i_a-1}$ is a normalisation factor and η_s, ia is the likelihood associated with the regime, s, which is given for each patch, i_a, in the Gaussian case by

$$\begin{array}{c}\hfill {\eta }_{s,i_a}={\displaystyle \sum _n^{{\theta }^2}\frac{1}{{\theta }^2}\frac{1}{\sqrt{2\pi }\sigma }exp[-\frac{{({\mathit{x}}_{i_a}^n-F_{i_a}\beta {\mathit{m}}^n-\mu )}^2}{2{\sigma }^2}],}\end{array}$$(3)

where θ gives the size in pixels of the planet model, m, and n is the pixel index within the patch.

Due to the dependence on the previous step, the RSM approach relies on an iterative algorithm to estimate the probability for every pixel contained in the selected annulus, a. Once all the probabilities of being in the planetary regime, ξ_1, ia, have been computed, those probabilities are averaged along the time axis to eventually provide the final RSM probability map. These different steps summarise the RSM approach, which was first proposed in Dahlqvist et al. (2020) and where a more detailed description of the algorithm can be found.

3. Using forward models in RSM

The original RSM map relies on an off-axis PSF to model the planetary signal. A promising development of the current method would be to take into account, via a forward model, the effects of the PSF subtraction techniques on the planetary signal. Indeed, most PSF subtraction techniques lead to distortions of the planetary signal, such as over-subtraction and self-subtraction (Pueyo 2016). Over-subtraction is attributable to quasi-static speckles inside the set of reference frames, while self-subtraction is due to the presence of the planetary signal itself inside the same set of reference frames. The signature of self-subtraction is specific to planetary candidates, as quasi-static speckles coming from the optical train do not rotate with the field. Because of field rotation, the evolution of the reference frames composition leads to the appearance of a negative wing travelling in time from one side of the planet to the other in the azimuthal direction. The temporal motion of this negative wing should therefore help disentangling a planetary candidate from a bright speckle.

We investigate, in this section, two forward model versions of the RSM map relying on the KLIP and LOCI algorithms. Both algorithms can accommodate an analytical estimation of the forward-modelled PSF, which is not the case for other ADI-based techniques such as NMF and LLSG. This avoid the complex task of choosing the fake companion intensity when constructing a forward-modelled PSF empirically by comparing an initial cube of residuals and one in which a fake companion has been injected.

3.1. KLIP-based forward modelling

Karhunen-Loéve image processing (KLIP) is a popular speckle subtraction technique first proposed by Soummer et al. (2012) and further improved by Pueyo (2016) who developed its forward model version. Similarly to PCA, KLIP estimates the reference PSF via a low-rank approximation of a reference library built to limit the impact of potential planetary signal on the speckle field estimation. For each frame of an ADI sequence, the KLIP algorithm computes the directions of maximal variance from the reference library. It keeps the principal components up to a rank, K, that is smaller than the dimension of the reference library, discarding the higher order modes that should contain more of the planetary signal. The principal components are found via a decomposition of the covariance matrix of the mean-subtracted reference frames R via the Karhunen-Loéve transform. They are given by:

$$\begin{array}{c}\hfill {\mathit{Z}}_K=\sqrt{{\mathbf{\Lambda }}^{-1}}{\mathit{V}}_K\mathit{R},\end{array}$$(4)

with Λ = diag(μ₁, μ₂, …, μ_k)^⊤ as the diagonal matrix with the eigenvalues of the image-to-image sample covariance matrix RR^⊤ with R the mean-subtracted reference library matrix and V_K = [v₁,v₂,…,v_k] as the respective eigenvectors up to the order K ≤ N_R, N_R being the number of images in the reference library used to compute the reference PSF (see Appendix A for a summary of all the variables used in the KLIP-based forward modelling). The reference PSF is then found by projecting the initial science image, i, or a subsection of this science image, onto the selected principal components. The reference PSF is subtracted from the science image yielding the residual image x as follows:

$$\begin{array}{c}\hfill \mathit{x}=\mathit{i}-{\mathit{Z}}_K^{\top }{\mathit{Z}}_K\mathit{i}.\end{array}$$(5)

The selection of the reference library is done via the definition of a minimal field of view (FOV) rotation between the science image, i, and the set of selected reference images, R. The minimal FOV rotation should be large enough to limit the distortion due to the planetary signal contained in the library (see Marois et al. 2010) but not too large so that it is possible to keep a sufficient correlation between the speckle field contained in the science image and the reference library. Pueyo (2016) proposed to model the distortion via an analytical expansion of the principal components to account for the presence of planetary signal inside the reference library. In the case of self-subtraction, the planetary signal appears in the principal components estimation via the covariance matrix and therefore the distortion is non linear as a result of the projection of m onto the perturbed components ΔZ_K. In contrast, the distortion due to over-subtraction is linear in m as it is defined as the projection of the planetary signal, m, on the unperturbed components, Z_K. The forward model of the planetary signal considers both type of subtractions as follows:

$$\begin{array}{c}\hfill \mathit{p}=\mathit{m}-{\mathit{Z}}_K^{\top }{\mathit{Z}}_K\mathit{m}-({\mathit{Z}}_K^{\top }\mathrm{\Delta }{\mathit{Z}}_K+{({\mathit{Z}}_K^{\top }\mathrm{\Delta }{\mathit{Z}}_K)}^{\top })\frac{\mathit{i}}{\beta },\end{array}$$(6)

where m represents the normalised planetary signal before reference PSF subtraction, typically the instrument off-axis PSF, and p is the forward model of the planet after subtraction. The second term on the right provides the over-subtraction of the point source while the third term gives the self-subtraction due to rotation via ΔZ_K (see Pueyo 2016, for the detailed derivation of this expression).

Having documented the estimation of the forward-modelled planetary signal and of the cube of residuals, we now consider how to include these elements in the RSM map framework. We rely on an annulus-wise estimation but in contrast with annular PCA or LLSG (Gomez Gonzalez et al. 2017), we do not estimate the speckle field for consecutive non-overlapping annuli. We estimate instead a specific speckle field for every radial distance, a. The self-subtraction wings appearing azimuthally, the brightest part of the planetary signal is contained in an annulus segment of width equal to one full width at half maximum (FWHM). The selected annulus with a width of one FWHM is centred on a and shifted by one pixel between each radial distance a instead of being shifted by one FWHM in the case of annular PCA or LLSG. This approach simplifies the forward model PSF estimation, provides a more accurate estimation of the speckle field, and avoids any non-linearities due to transitions between annuli.

The reference PSF and the forward-modelled PSF computation is done via Eqs. (5) and (6), respectively. The resulting forward-modelled PSF is then derotated and cropped to form a set of patches, p_ia, where all the elements outside the selected annulus segment is set to zero as can be seen in Fig. 1. Several crop sizes, from one to two FWHM, are tested in the next section. Once injected in the expression of the likelihood given in Eq. (3), this allows us to focus the model only on the region where strong intensity variations occur. The expression of the likelihood of being in the state, s, for every patch, i_a, becomes (in the Gaussian case):

$$\begin{array}{c}\hfill {\eta }_{s,i_a}={\displaystyle \sum _n^{{\theta }^2}\frac{1}{{\theta }^2}\frac{1}{\sqrt{2\pi }\sigma }exp[-\frac{{({\mathit{x}}_{i_a}^n-F_{i_a}\beta {\mathit{p}}_{i_a}^n-\mu )}^2}{2{\sigma }^2}],}\end{array}$$(7)

Fig. 1. KLIP forward-modelled PSF for a NACO β Pictoris ADI sequence taken at a distance of 2λ/D, cropped at two FWHM and summed along the time axis.

where x_ia are the derotated and cropped patches obtained from the residual image, x.

When relying on a forward-modelled PSF, the intensity parameter β directly provides an estimation of the planet luminosity, which may be helpful for characterizing the planetary candidate beyond its detection. Two methods are considered for estimating this intensity parameter. The first method is similar to the one proposed in Dahlqvist et al. (2020), with β defined as a multiple of the estimated standard deviation of the pixel intensity in the annulus, β = δσ. In this case, the standard deviation σ is estimated empirically by considering all the frames and an annulus with a width of one FWHM centred on the annulus of interest, a. The δ parameter is defined via the maximisation of the total likelihood of the annulus. The second method relies on the definition of the intensity via a Gaussian maximum likelihood (see Cantalloube et al. 2015; Ruffio et al. 2017, for more details) before the computation of the RSM map itself, which allows the use of the following analytical form for the flux parameter β:

$$\begin{array}{c}\hfill \stackrel{\sim }{\beta }=\frac{{\sum }_j^T{\mathit{i}}_j^{\top }{\mathit{p}}_j/{\sigma }_j}{{\sum }_j^T{\mathit{p}}_j^{\top }{\mathit{p}}_j/{\sigma }_j},\end{array}$$(8)

with the standard deviation σ_j computed separately for each frame by considering an annulus with a width equal to one FWHM centred on a.

The main advantage of this second method is the simplicity of the intensity computation. It provides also a specific intensity for each pixel, which may help the algorithm to differentiate bright speckles from planetary candidates. The only drawback of this approach is that it makes the assumption of a Gaussian residuals distribution, which is not always the case especially near the host star. As for the parametrisation of the RSM likelihood function in this second case, the standard deviation used for the flux estimation is taken and the computation of the mean is also done frame-wise using the same procedure. Both approaches are investigated in Sect. 4.

3.2. LOCI-based forward modelling

The second PSF subtraction technique that we propose for investigating the RSM framework is a forward model version of LOCI, the locally optimised combination of images (Lafreniere et al. 2007). Specifically, LOCI relies on a linear combination of reference images to model the speckle field in a given science image, i. As for the KLIP algorithm, the definition of the reference library is based on the definition of a minimal FOV rotation between the frames composing the reference library and the selected science image, that is, a minimal distance by which a potential point source in the science image would be displaced in the frames composing the reference library.

Besides this angular distance, LOCI relies on the definition of two different subsections within the ADI sequence. A first section, O_K, is used for the computation of the linear combination factors, while a second smaller subsection R_K is selected for the speckle field subtraction². The use of a larger section O_K aims to reduce the weight of the potential planetary candidate in the estimation of the linear combination. Once the reference library is defined, the computation of the linear combination factors is simply done via the minimisation of the sum of squared residuals (Lafreniere et al. 2007):

$$\begin{array}{c}\hfill {ϵ}^2={\displaystyle \sum _{j=0}^{N_p}{[{\mathit{o}}_i^j-\sum _k^{N_R}{c}_k{\mathit{o}}_k^j]}^2,}\end{array}$$(9)

with o_i the section of the frame for which a model of the speckle field is computed via the factors c_k, and o_k the section of the reference frame k (see Appendix A for a summary of all the variables used in the LOCI-based forward modelling). The minimum of this last expression has an analytical form obtained by setting all the partial derivatives with respect to c_k equal to zero, which is equivalent to solving a simple system of linear equations of the form Ax = b given by:

$$\begin{array}{c}\hfill {\displaystyle \sum _k^{N_R}{c}_k\left[\sum _{j=0}^{N_p}{\mathit{o}}_l^j{\mathit{o}}_k^j\right]=\left[\sum _{j=0}^{N_p}{\mathit{o}}_l^j{\mathit{o}}_i^j\right],}\end{array}$$(10)

which holds ∀l ∈ K.

Once obtained, the factors, c_k, are multiplied by the subsection of the reference frames, r_k, and subtracted from the subsection of the science image, i, to get the residual, x, as follows:

$$\begin{array}{c}\hfill \mathit{x}=\mathit{i}-{\displaystyle \sum _k^{N_R}{c}_k{\mathit{r}}_k,}\end{array}$$(11)

i defined here as the subsection of the science image corresponding to the reference frames, r_k.

The forward model of the planetary signal is easily computed using the same factors, c_k, and the planetary signal, m:

$$\begin{array}{c}\hfill \mathit{p}={\mathit{m}}_i-{\displaystyle \sum _k^{N_R}{c}_k{\mathit{m}}_k.}\end{array}$$(12)

As in the case of the KLIP forward model, for both the residual image and the forward model PSF, an annulus with a width of one FWHM is used to focus on the region where the planetary signal is the most visible. The forward model PSF, p, and the residuals images, x, are again derotated and cropped to from the time series, x_ia and p_ia. The two methods used to estimate the flux parameter β are again considered.

3.3. Forward model RSM map: summary

We briefly summarise the main steps for the RSM map estimation when relying on the LOCI and KLIP forward-model approach as follows:

1. Compute the residuals for an annulus centred on a for each frame using the KLIP or LOCI procedure.

2. Compute the PSF forward model for every frame and every position within the annulus, a.

3. Derotate the resulting annuli and crop the forward-modelled PSF and science image to form the times series, x_ia and p_ia.

4. Estimate the mean and variance of the residuals for every frame, considering the annulus of width equal to one FWHM centred on a.

5. Using the iterative procedure described in Sect. 2, estimate, ξ_1, ia for each index, i_a, using the forward model version of the likelihood (see Eq. (8)).

6. Repeat steps 1 to 5 for every annulus.

7. Average the resulting probability matrix along the time axis to obtain the final RSM detection map.

4. Performance assessment of a forward-modelled RSM map

4.1. Data sets

We propose a reliance on data sets provided by three different instruments to assess the performance of the two forward model versions of the RSM map. The two first ADI sequences were acquired with two instruments of the Very Large Telescope (VLT), NACO, and SPHERE, while the third sequence was acquired with the LMIRCam instrument of the Large Binocular Telescope (LBT). This choice of data sets aims to investigate the behaviour of the algorithm when facing different noise profiles generated by a variety of instruments.

The first data set is an ADI sequence on β Pictoris and its planetary companion β Pictoris b obtained in L band in January 2013 with NACO in its AGPM coronagraphic mode (Absil et al. 2013). The ADI sequence is composed of 612 individual frames obtained by integrating 40 successive individual exposures of 200 ms. Every third frame was selected here to reduce the computation time, resulting in a final cube of 204 frames. The parallactic angle ranges from −15° to +68°.

The second ADI sequence focuses on 51 Eridani. It was obtained in K1 band in September 2015 with the SPHERE-IRDIS instrument, using an apodised pupil Lyot coronagraph (Samland et al. 2017). The data set regroups 194 frames pre-processed using the SPHERE Data Center pipeline (for more details about the reduction see Delorme et al. 2017; Maire et al. 2019). The integration time is 16 s and the parallactic angle ranges from 297° to 339°.

The last data set is an ADI sequence on HD183324 produced by the LMIRCam instrument of the LBT. The images were obtained in October 2018 in L′ band without coronograph using a single telescope. The pre-processed data set³contained 1394 frames with integration time of 109 ms. They were binned over 10 successive individual exposures to reduce the computation time, leading to 139 frames with an integration time of 1.09 s. The parallactic angles range from −13° to −39°. A region with a radius of one arcsecond is considered for all three data sets, which corresponds to around 16 λ/D for the SPHERE data set and 8 λ/D for the NACO and LMIRCam data sets.

4.2. Results

The performance assessment of the two forward model versions of the RSM map is done via the estimation of ROC curves. In contrast to the ROC curve usually used for assessing the performance of binary classifiers, the false positive rate (FPR) is replaced by the number of false positive (FP) for the entire frame, averaged over the number of test data sets used for the ROC curve computation (see Dahlqvist et al. 2020 and Gomez Gonzalez et al. 2018). Synthetic data sets are generated based on the three selected ADI sequences by injecting fake companions at two different angular separations to account for the radial evolution of the noise profile. The known companions and some bright disk structures for the β Pictoris data set were removed via the negative fake companion technique (Lagrange et al. 2010) prior to generating the synthetic data sets. The fake companions, which are simply defined as the normalised off-axis PSF, are injected at 16 different position angles with five different flux values for a given angular separation. This allows us to test the sensitivity of the forward model RSM map to different contrasts and mitigates the impact of local speckles on the estimation of the ROC curves. The contrasts used for the three ADI sequences are given in Table 1. The relatively low contrasts used for the LMIRCam data set arise from the short integration time, the low number of frames after binning as well as the small angular rotation, all of which affect the performance of the PSF subtraction techniques. This provides an interesting way of exploring the algorithm performance in different HCI regimes.

We consider a true positive (TP) for a given threshold to be a peak value above the threshold in a circle with a diameter of one FWHM centred on the position of the injected fake companion. A value above the selected threshold at any other location is considered as a FP. In order to avoid double counting, we impose the condition that peak values outside the fake companion region should be separated by a minimal distance of one FWHM.

4.2.1. KLIP-FM RSM map

The forward model version of the KLIP algorithm was developed along a Gaussian matched filter to detect potential planetary candidates in the cube of residuals using the PSF forward model. We propose therefore to compare the performance of the forward model RSM map with the performance of the KLIP forward model matched filter (KLIP-FMMF) developed by Ruffio et al. (2017). We include additionally the original RSM map applied on the cube of residuals generated by the KLIP PSF subtraction techniques and the S/N map obtained with KLIP. In the last case, the S/N map is generated annulus-wise using the procedure of Mawet et al. (2014), which estimates the S/N by comparing the flux inside a one FWHM aperture centred on the considered pixel with the flux of all the other apertures included in the annulus. This procedure implements a small-sample statistics correction, relying on a student t-test to determine the S/N. Increasing S/N or probability thresholds are applied to generate the different ROC curves for all the considered methods.

The parameters of the KLIP algorithm, namely, the number of principal components and minimum FOV rotation, were selected to optimise the ROC curves⁴ for the two considered angular separations. A single set of parameters was defined for each data set. The number of principal components was set to 20 for the SPHERE and NACO data sets while a value of 18 principal components was chosen for the LMIRCam data set. The FOV rotations expressed in terms of FWHM are respectively 0.5, 0.3, and 0.3.

We start by considering for all three data sets the impact of the selected crop size for the forward-modelled PSF used in the KLIP-FM RSM map. As can be seen from Fig. 2, the larger crop sizes seem to outperform the crop size of one FWHM for the small angular separation while the reverse is true for the large angular separation. However we observe a much larger gap for the largest separation, especially in the case of the NACO data set. This may be explained by the reduced self-subtraction observed at large angular separations, the movement of potential astrophysical signals increasing linearly with the angular separation. This implies that the brightness of the negative lobes appearing on the sides of the main peak reduces with the angular separation. This makes the larger crop sizes unnecessary and more prone to speckle noise.

Fig. 2. ROC curves for the NACO, SPHERE, and LMIRCam data sets, with the KLIP-FM RSM map using respectively a crop size for the froward modelled PSF of 5 (red), 7 (blue), 9 (green), 11 (orange) pixels (FWHM ≈ 5 pixels for all three data sets). (a) NACO at 2λ/D, (b) SPHERE at 2λ/D, (c) LMIRCam at 2λ/D, (d) NACO at 8λ/D, (e) SPHERE at 8λ/D, (f) LMIRCam at 8λ/D.

Considering the previous results, we select for the NACO data set a crop size of five pixels (one FWHM) and a crop size of seven pixels for the other two data sets. This provides a good performance trade-off between small and large angular separations. As can be seen from Fig. 3, a crop size of one FWHM still captures part of the negative wing azimuthal translation. The results for the two versions of KLIP-FM RSM map, as well as KLIP-FMMF, KLIP RSM map and the KLIP using S/N map are given in Fig. 4. We see from these plots that, at small separation, the KLIP-FM RSM map seems to slightly outperform KLIP-FMMF, while the reverse is true at large separation. The KLIP approach using S/N map has a higher ability to detect faint companions at large radial distances but it is no match to the other methods at small separations. The KLIP RSM map provides surprisingly good results, being often the closest to the KLIP S/N map for large separations and being relatively close to KLIP-FM RSM and KLIP-FMMF at 2λ/D from the host star. A combination of both KLIP RSM and KLIP-FM RSM could be interesting to keep the high sensitivity of KLIP-FM RSM at close separations while improving the sensitivity at larger radial distances. It seems also clear from Fig. 4 that the version of the KLIP-FM RSM map using the Gaussian approximation for estimating the flux parameter β (Eq. (8)) outperforms the one relying on the maximum likelihood approach used in Dahlqvist et al. (2020), providing in all cases equivalent or better results. We also tested the maximum likelihood based approach and the Gaussian approximation with the KLIP RSM map with similar results (see Appendix B for a comparison between the two approaches in the case of KLIP RSM), demonstrating the efficiency of this new way of estimating β on top of its faster estimation.

Fig. 3. KLIP forward-modelled PSF for the NACO β Pictoris ADI sequence taken at a distance of 2λ/D with the same azimuthal orientation as in Fig. 1 and cropped at one FWHM. The two images correspond respectively to the (a) first and the (b) last frame.

4.2.2. LOCI-FM RSM map

Turning to the RSM LOCI FM map, the tolerance level for the square-residuals minimisation and the minimum FOV rotation were also selected to provide the best overall performance. A tolerance of 9 × 10⁻³ was chosen for NACO and SPHERE and a tolerance of 1 × 10⁻² for LMIRCam. The minimum FOV rotations are respectively, 0.6, 0.2 and 0.2 FWHM. The analysis of the ROC curves obtained with different crop sizes leads to similar conclusions to the case of KLIP-FM RSM. The crop size of one FWHM performs better, in a global sense, even though larger crop sizes lead to slightly better results at small angular separations. The ROC curves corresponding to the crop sizes performance comparison are presented in Appendix C. Regarding the performance of LOCI-FM RSM, the results in Fig. 5 demonstrate again the interest of the Gaussian maximum likelihood to define the flux parameter β mainly for the largest separation. Both LOCI-FM RSM and LOCI RSM outperform clearly the LOCI S/N map for the 2λ/D angular separation, while the reverse is true for the 8λ/D angular separation. The ordering is similar to that of the KLIP case with LOCI-FM RSM leading at 2λ/D and LOCI RSM being closer to the LOCI S/N map at 8λ/D, which again seems to favour a combination of both the LOCI and LOCI-FM to benefit from their respective strength. The search for an optimal mix between the different PSF subtraction techniques is investigated in the next section.

5. Optimal PSF subtraction techniques selection

Having demonstrated the added value of the forward model versions of the RSM map, at least at small angular separations, we are now left with five different PSF subtraction techniques (annular PCA, KLIP, NMF, LLSG, LOCI) plus two forward model versions to generate the RSM maps. Given that the Annular PCA and KLIP are relatively close in their definition, we decided to focus solely on KLIP, as preliminary results demonstrated their similarities in terms of performance and their non-complementarity. We address in this section the difficult question of optimally selecting these PSF subtraction techniques to optimise the overall performance of the resulting RSM maps. In particular, we investigate the dependence of the optimal combination on the instrument and radial distance. We rely again on ROC curves to assess the performance of the various combinations we considered.

In order to speed up the multiple RSM map estimations, we slightly modified the original RSM map procedure as presented in Dahlqvist et al. (2020), with, however, no impact on the final outcome of the algorithm. We divided the procedure into two separate steps, the first one being the estimation of the likelihood provided in Eq. (3) and the second one the estimation of the probability of being in the planetary regime given by Eq. (2). A separate likelihood cube is estimated for every considered PSF subtraction technique for the entire set of annuli. Some of these likelihood cubes are then stacked along the time axis depending on the selected combination. The probabilities are eventually estimated annulus-wise for every pixel of every frame and averaged along the time axis to generate the final probability map⁵. This allows us to estimate only once the likelihood cubes for the different PSF subtraction techniques, with the second step, which is also the fastest, as the only one to be repeated for each combination. The parametrisation of the underlying PSF subtraction techniques were selected to maximise the overall performance for each data set.

The ROC curve computation follows the same procedure as in the previous section with a region of one arcsecond considered for all three data sets and the curves being computed for the two same angular separations. The main characteristics of the ROC curves for the 16 selected combinations of PSF subtraction techniques may be found in Appendix D. The two parameters we used to select the best combinations, are the maximum TPR reached without any FP, and the average number of FPs inside the entire frame at TPR = 1⁶. The first parameter is the most important one, as it gives clues about the highest contrast the algorithm can reach without any false detection. The second gives a measure of the number of bright background structures that have not been properly treated by both the PSF subtraction techniques and the RSM algorithm.

The results presented in Appendix D⁷ show large differences in terms of performance between the considered combinations highlighting the importance of the PSF subtraction techniques selection. Based on the two metrics introduced in the previous paragraph, we selected the four best combinations for each data set. The ROC curves for the two considered angular separations are presented in Fig. 6. When comparing these curves with the ones in Figs. 4 and 5, we see that the improvement of the RSM map performance occurs mainly at larger separations when considering multiple PSF subtraction techniques. The ROC curves are indeed very close to the ones obtained with the KLIP-FM RSM and LOCI-FM RSM for the 2λ/D radial distance, while the gap is much wider for the 8λ/D. Apart from the combinations 7 and 14 (see Fig. D.1), which are selected for multiple data sets, the other combinations are specific to each data set. This seems to demonstrate that the selection of an optimal combination should be done at least on an instrument-specific basis. The definition of a single optimal combination for the entire set of annuli seems also difficult, as we often observe that higher performance at short separations goes hand in hand with lower performance at large separations. A last element to consider for the selection of the optimal combination is the threshold value for which the first false positive is observed, which should be as small as possible since large values imply the presence of bright structures in the probability map. We conclude from Fig. 6 that in a global sense, the best combinations are the combinations 3, 11, and 10 for the NACO, SPHERE, and LMIRCam data sets, respectively. Looking at these three combinations, we see that they share a common structure, being composed of the LLSG and NMF PSF subtraction techniques, in addition to a LOCI- or KLIP-based PSF subtraction technique. The performance of this particular combination is probably due to the differences between these PSF subtraction techniques in terms of residuals noise profile. These differences should help to better average out the speckle noise via the RSM algorithm. This structure therefore appears to be an interesting starting point when studying a new data set. We further characterise these three combinations in Sect. 7 by estimating their contrast curve.

Fig. 4. ROC curves for the NACO, SPHERE and LMIRCam data sets, with respectively, the KLIP-FM RSM map using the Gausssian maximum likelihood for the pre-optimisation of the flux parameter β (red), the KLIP-FM RSM map with no flux pre-optimisation (NF), which relies on the standard maximum likelihood used in the original RSM map for the estimation of flux parameter β (blue), the forward model matched filter KLIP-FMMF (green), the RSM map using KLIP (orange) and KLIP using the standard S/N map (black). (a) NACO at 2λ/D, (b) SPHERE at 2λ/D, (c) LMIRCam at 2λ/D, (d) NACO at 8λ/D, (e) SPHERE at 8λ/D, (f) LMIRCam at 8λ/D.

Fig. 5. ROC curves for the NACO, SPHERE and LMIRCam data sets, with, respectively, the LOCI-FM RSM using the Gausssian maximum likelihood for the pre-optimisation of the flux parameter β (red), the LOCI-FM RSM map with no flux pre-optimisation (NF), which relies on the standard maximum likelihood used in the original RSM map for the estimation of flux parameter β (blue), the RSM map using LOCI (orange) and LOCI using the standard S/N map (black). (a) NACO at 2λ/D, (b) SPHERE at 2λ/D, (c) LMIRCam at 2λ/D, (d) NACO at 8λ/D, (e) SPHERE at 8λ/D, (f) LMIRCam at 8λ/D.

Fig. 6. ROC curves for the NACO, SPHERE and LMIRCam data sets, with the four best combinations of PSF subtraction techniques used to generate the RSM map algorithm. (a) NACO at 2λ/D, (b) SPHERE at 2λ/D, (c) LMIRCam at 2λ/D, (d) NACO at 8λ/D, (e) SPHERE at 8λ/D, (f) LMIRCam at 8λ/D.

The results presented in this section demonstrate the dependence of the optimal set of PSF subtraction techniques on the instrument providing the ADI sequence, but also on the angular separation, although a common underlying structure could be seen. A larger set of ADI sequences would be needed to determine whether a single optimal set of PSF subtraction techniques may be identified for a given instrument or even multiple instruments, which could be very helpful when dealing with large surveys.

6. Forward-backward model

In this section, we discuss an additional improvement of the original RSM map by considering a forward-backward approach for the estimation of the probability, ξ_1, ia. The current approach relies solely on past observations to construct the cube of probabilities while the entire cube of residuals is available for the estimation, that is, of both past and future observations. We propose therefore to replace the current forward approach by a forward-backward approach, which considers both past and future observations. This method computes two separate sets of probabilities, the forward probabilities as done in the original RSM framework:

$$\begin{array}{c}\hfill {\xi }_{1,i_a}^f={\displaystyle \sum _{q=0}^1\frac{{\eta }_{1,i_a}{p}_{q,1}{\xi }_{q,i_a-1}^f}{{\sum }_{q=0}^1{\sum }_{s=0}^1{\eta }_{s,i_a}{p}_{q,s}{\xi }_{q,i_a-1}^f},}\end{array}$$(13)

but also the backward probabilities, which rely on the probability estimated at index i_a + 1 instead of index i_a − 1 to compute the probability at current index i_a as:

$$\begin{array}{c}\hfill {\xi }_{1,i_a}^b={\displaystyle \sum _{q=0}^1\frac{{\eta }_{1,i_a}{p}_{q,1}{\xi }_{q,i_a+1}^b}{{\sum }_{q=0}^1{\sum }_{s=0}^1{\eta }_{s,i_a}{p}_{q,s}{\xi }_{q,i_a+1}^b}.}\end{array}$$(14)

Once both sets have been estimated, the final probabilities are obtained by multiplying the two sequences of probabilities. A normalisation factor is applied, making sure that the total probability equals 1 for every index i_a. The final probabilities are therefore given by:

$$\begin{array}{c}\hfill {\xi }_{1,i_a}=\frac{{\xi }_{1,i_a}^f{\xi }_{1,i_a}^b}{{\sum }_{s=0}^1{\xi }_{s,i_a}^f{\xi }_{s,i_a}^b}\end{array}$$(15)

Because the RSM features a short-term memory, the probability of being in the planetary regime builds up when we get closer to the planetary signal but with a small latency. As can be seen from Fig. 7, this latency leads to a shift of the main peak towards the future for the forward approach and towards the past for the backward approach. When relying on the forward-backward approach these shifts cancel out and the main peak is centred on the true position of the planetary signal. The forward-backward approach should therefore allow to reach a higher precision in terms of exoplanet astrometry.

Fig. 7. Evolution of the probabilities for the forward, backward, and forward-backward approaches using KLIP, around the location of a planetary candidate injected in the 51 Eridani data set (radial distance of 4λ/D with a contrast of 3.76 × 10−5).

In order to investigate the ability of both approaches to derive an accurate astrometric measurement for the detected planetary signal, we propose performing a series of simulations based on the SPHERE data set. We study the evolution of the astrometric precision for a range of contrasts, considering again a radial distance of 2λ/D and 8λ/D. As done in the previous sections, we base our simulation on synthetic data sets, on the basis of which we apply a KLIP RSM map using the forward and forward-backward version of the RSM algorithm. The negative fake companion (NEGFC) method (Lagrange et al. 2010; Marois et al. 2010; Wertz et al. 2017) is also applied on the synthetic data sets, allowing for a comparison with a technique dedicated to astrometry⁸. For each radial distance, we inject fake companions at 16 different position angles. The set of considered contrasts are computed based on the KLIP RSM contrast curve, estimated using the approach proposed in the next section (see Fig. E.1). We define two sets of contrasts ranging from one to six times the sensitivity limit at the considered radial distance in Fig. E.1, with a step size of 0.5. In the case of the RSM map, the computation of the position is done by fitting a two-dimensional Gaussian to the detected planetary signal. The astrometric error bars for the three considered methods are computed as the root mean squared (rms) position error between the obtained position and the injected fake companion true position, averaged over the two axes. The rms is estimated over the 16 fake companions injected at each radial distance, for every contrast.

The results from Fig. 8 demonstrate clearly the ability of the forward-backward approach to decrease the position error compared to the original forward approach. As can be seen from Fig. 8, the RSM forward-backward approach performs better than the NEGFC approach at large radial distances and for high contrasts. However, for lower contrast, the RSM forward-backward approach reaches a noise floor around 4 mas, higher than the noise floor of the NEGFC approach, which is between 1.5 and 2 mas. This higher noise floor may be explained partly by the profile of the planetary signal in the RSM map. As can be seen from Eq. (2), the RSM approach response to a planetary signal is non linear and dependent on neighbouring pixels, leading potentially to asymmetries in the azimuthal direction, even in the forward-backward case. The algorithm architecture also leads to non-linearities along the radial axis because of the annulus-wise probabilities computation. Finally, as can be seen from Fig. 7, the forward-backward approach reduces the amplitude of the planetary signal within the probability map. All these elements affect the Gaussian fit and therefore the astrometric precision that the RSM algorithm can reach. Nevertheless, as demonstrated by the results from Fig. 8, the RSM forward-backward approach can reach a higher astrometric precision, especially at large radial distances and high contrasts. This is due to the better ability of the RSM algorithm to detect faint companions. It is also worth noting that the computation time is much lower when using the RSM map than with the NEGFC approach or the more advanced Markov Chain Monte Carlo version of NEGFC approach. The RSM forward-backward approach provides therefore a good first estimate, especially for high contrast targets, which can then be used, for lower contrasts, as an initial position for more advanced astrometry techniques.

Fig. 8. RMS position errors averaged over the two axes expressed in mas for the KLIP RSM map using the forward (blue) and forward-backward (red) versions of the RSM approach and for the NEGFC approach (orange). The two graphs show the dependence of the averaged RMS position error on the contrast for respectively a radial distance of (a) 2λ/D and (b) 8λ/D.

Another advantage of the forward-backward approach lies in its ability to reduce the background speckle noise and smooth the probability curve, the noise being treated differently by the forward and backward components. Looking at Fig. 9, we see that the level of the residual speckle noise has reduced drastically for the LMIRCam data set, the brightest speckle probability decreasing by around 40 percent. However this noise reduction comes along with a reduced brightness of the planetary signal. This reduced brightness is also illustrated in Fig. 7, where the peak value obtained with the forward-backward approach is lower than the one obtained with the two other approaches, leading to a reduced detection threshold. This reduction of the planet signal strength does not impact the performance of the forward-backward approach in terms of ROC curve, with similar results for both the forward and forward-backward approaches. The forward-backward approach outperforms even sightly the forward approach at small separation (see Appendix B for a comparison between the original RSM approach and the forward-backward version in the case of KLIP).

Fig. 9. RSM map generated with forward and forward-backward approach for the LMIRCam data set using KLIP with 18 principal components, a FOV rotation expressed in terms of FWHM of 0.3 and a Gaussian distribution. A square root-based scale has been used to increase the background speckle noise brightness. (a) Original RSM map and (b) Forward-Backward RSM map.

7. Contrast curve

This last section is devoted to the estimation of contrast curves based on our RSM framework. When relying on probability map for exoplanet detection, we cannot use the traditional procedure to compute the contrast curve. This procedure selects the contrast corresponding to a TPR of 0.5 for a given probability of observing a FP (Jensen-Clem et al. 2017). A 5σ threshold is usually chosen, which corresponds to a 3 × 10⁻⁷ false alarm probability under a Gaussian noise hypothesis. In the case of the RSM map, the probability provided as an output by the algorithm is a non linear function of the underlying likelihoods, past observations, and the transfer matrix (see Eq. (2)), which precludes us from defining a similar probability threshold for the FP. We therefore rely on a procedure similar to the one we used for the ROC curve computation. Considering a given flux, we inject fake companions at different position angles for a given radial distance. The resulting synthetic data sets allow for the computation of a threshold-dependent TPR. The next step is to select a convenient threshold. As it is not possible to reach a 5σ confidence empirically, we select the threshold that leads to the first detection of a FP within the entire frame (as done in Jensen-Clem et al. 2017, but considering the entire frame instead of the selected annulus). This threshold definition provides a direct link with the background noise level while avoiding the usual shortcomings of standard contrast curves, namely, the noise Gaussianity assumption and the definition of the region used to empirically estimate the first two moments of the noise.

We repeat the steps presented in the previous paragraphs on sets of fake companions injected with different flux values until a TPR of 0.5 is found when detecting the first FP. We use an iterative procedure based on linear interpolations to minimise the number of flux values to be considered before reaching the TPR of 0.5. The procedure starts with the estimation of the TPR for a selected pair of upper and lower flux values (step 1 in Fig. 10, respectively below and above the TPR of 0.5 preferably). A linear interpolation then allows the determination of a best guess for the next flux value to be estimated (step 2). Once the TPR is estimated for this new flux value, a new linear interpolation is performed with this last flux value as upper bound (resp. lower bound) if the TPR is above 0.5 (resp. below 0.5) keeping the previous lower flux value (resp. upper flux value) (step 3). The procedure is repeated until a TPR ∈[0.45,0.55] is found using a tolerance interval with a size of 0.1. These steps are summarised in Fig. 10.

Fig. 10. Procedure for the estimation of the contrast curve, the orange star corresponds to the linear interpolation, the black star to the corresponding TPR, the blue stars to the two previous estimations used for the linear interpolation and the green star to the flux corresponding to a TPR of 0.5. The red curve makes the link between the flux and the true positive rate in the case of a single false positive while the green line shows the tolerance’s interval of [0.45,0.55].

We estimated contrast curves for the three data sets, considering the original RSM map using the optimal combination found in Sect. 7, the forward-backward version of the algorithm, and the simple KLIP S/N map. We computed the contrast from a radial distance of 2λ/D up to 8λ/D using a step-size of one λ/D and then up to 16λ/D for the SPHERE data set with a step-size of two λ/D. The resulting contrast curves may be found in Fig. 11. These contrast curves confirm the findings made with the ROC curves, the RSM map reaching in most cases a significantly higher contrast for small angular separations and providing comparable contrast at larger separations. The gap between the contrast curves of the KLIP S/N map and the RSM map is significantly larger for the LMIRCam data set, highlighting the interest of the RSM approach for this instrument. However, the lower number of frames and the smaller angular rotation for this particular data set may explain these results. A larger number of LMIRCam data sets is therefore needed to confirm these findings. We also find that the RSM map and the forward-backward version have a similar performance for all angular separations. Considering its higher precision and its lower background noise level, the forward-backward version seems to be a promising alternative to the standard RSM map.

Fig. 11. Contrast curves for the three data sets using the optimal combination of PSF subtraction techniques obtained in Sect. 5 (resp. 3, 11, and 10) with the region [2λ/D,16λ/D] considered for the SPHERE data set and the region [2λ/D,8λ/D] considered for the other two data sets to get the contrast in the first arcsecond. (a) NACO, (b) SPHERE, (c) LMIRCam.

8. Conclusion

In order to improve the RSM map sensitivity to faint exoplanets, we implemented several improvements at different levels. We first considered a forward-model version of the RSM map for two different PSF subtraction techniques, KLIP and LOCI. The computation of forward-modelled PSFs allowed to take into account PSF distortions due to the speckle field subtraction. We demonstrated the interest of the forward model at small angular separations through the estimation of ROC curves. These ROC curves highlight the radial dependence of the optimal PSF crop size, with larger crop sizes leading to better results for small angular separations while the reverse is true for smaller crop sizes. This could be explained by the fact that larger PSF crop sizes better account for the self-subtraction patterns, whose intensity decreases with radial distance as the relative movement of astrophysical signals increases with the distance to the host star.

In a second step, we investigated the question of the optimal selection of the PSF subtraction techniques to be included in the estimation of the RSM map. Relying again on ROC curves, we demonstrated the importance of this selection with sometimes large differences between the performance of the considered combinations. The optimal combinations seems to depend on the instrument, but also on the angular separation.

We proposed then an improvement directly related to the way probabilities are estimated within the RSM algorithm by replacing the original forward approach by a forward-backward approach. This approach relies on both past and future observations to compute a final probability. This allows us to deal more efficiently with background noise as speckles are not treated in the same way when relying on a forward and backward approach. Another advantage of the forward-backward version of the RSM map is its ability to extract more precisely the planetary astrometry.

We finally implemented a new framework to compute contrast curves in the case of probability maps, which cannot be treated as S/N maps. We kept the TPR of 0.5 while removing the Gaussian based noise threshold definition, replacing it with a threshold based on the detection of the first FP. Using a procedure based on successive linear interpolations, we were able to compute the contrast curve for the original RSM map and the forward-backward version using the optimal combination obtained previously and compare their performance with the one of a simple KLIP S/N map. The results highlighted the ability of the RSM map to detect fainter companions at small angular separations (below 400 mas) and the overall similar results for larger separations. The contrast curves also confirmed the interest of the forward-backward approach as it provides similar contrast curves while reducing the speckle background noise and giving a higher precision in terms of astrometry.

As the selection of the optimal combination of PSF subtraction techniques to be included in the RSM map as well as the parametrisation of these PSF subtraction techniques are relatively time-consuming, an interesting future development of the RSM map would be to create an automated parametrisation procedure that depends on angular separation.

Acknowledgments

This work was supported by the Fonds de la Recherche Scientifique - FNRS under Grant n° F.4504.18 and by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement n° 819155). We thank our colleagues A.-L. Maire for sharing the SPHERE data set and Arianna Musso-Barcucci for sharing the LMIRCam data set.

References

Appendix A: Mathematical notations for RSM map, KLIP FM and LOCI FM

Appendix B: Estimation of the β parameter and forward-backward model in the case of RSM KLIP

Fig. B.1. ROC curves for the NACO, SPHERE, and LMIRCam data sets, with respectively the KLIP RSM map using the Gaussian maximum likelihood for the pre-optimisation of the flux parameter β (red), the KLIP RSM map with no flux pre-optimisation (NF), which relies on the standard maximum likelihood used in the original RSM map for the estimation of flux parameter β (blue), and the KLIP RSM map using the forward-backward approach for the probability estimation instead of the original forward approach (green). (a) NACO at 2λ/D, (b) SPHERE at 2λ/D, (c) LMIRCam at 2λ/D, (d) NACO at 8λ/D, (e) SPHERE at 8λ/D, (f) LMIRCam at 8λ/D.

Appendix C: Crop size comparison for FM LOCI RSM

Fig. C.1. ROC curves for the NACO, SPHERE, and LMIRCam data sets, with the LOCI-FM RSM map using respectively a crop size for the froward modelled PSF of 5 (red), 7 (blue), 9 (green), 11 (orange) pixels (FWHM ≈ 5 pixels for all three data sets). (a) NACO at 2λ/D, (b) SPHERE at 2λ/D, (c) LMIRCam at 2λ/D, (d) NACO at 8λ/D, (e) SPHERE at 8λ/D, (f) LMIRCam at 8λ/D.

Appendix D: Performance comparison of PSF subtraction techniques combinations

Fig. D.1. Cumulative maximum TPR reached without any FP (blue), and cumulative average number of FPs inside the entire frame at TPR = 1 (red) for 16 different combinations of PSF subtraction techniques used to generate the RSM map and two radial distances. The dark-coloured bars are obtained at a radial distance of 2λ/D while the light-coloured bars correspond to 8λ/D for, respectively, the NACO (top), SPHERE (middle), and LMIRCam (bottom) data sets. The asterisks for some values of FP at TPR = 1 indicate that a TPR of 1 has not been reached at a distance of 8λ/D and that the smallest probability threshold was chosen instead (highest TPR). A high performance for a combination of PSF subtraction techniques corresponds to a tall blue histogram alongside a short red histogram.

Appendix E: Contrast curve for the SPHERE data set using KLIP RSM

Fig. E.1. Contrast curves for the SPHERE data set using KLIP RSM with 20 principal components and a FOV rotation expressed in terms of FWHM of 0.3. The region [2λ/D,16λ/D] has been considered to get the contrast in the first arcsecond.

All Tables

All Figures

	Fig. 1. KLIP forward-modelled PSF for a NACO β Pictoris ADI sequence taken at a distance of 2λ/D, cropped at two FWHM and summed along the time axis.
In the text

	Fig. 2. ROC curves for the NACO, SPHERE, and LMIRCam data sets, with the KLIP-FM RSM map using respectively a crop size for the froward modelled PSF of 5 (red), 7 (blue), 9 (green), 11 (orange) pixels (FWHM ≈ 5 pixels for all three data sets). (a) NACO at 2λ/D, (b) SPHERE at 2λ/D, (c) LMIRCam at 2λ/D, (d) NACO at 8λ/D, (e) SPHERE at 8λ/D, (f) LMIRCam at 8λ/D.
In the text

	Fig. 3. KLIP forward-modelled PSF for the NACO β Pictoris ADI sequence taken at a distance of 2λ/D with the same azimuthal orientation as in Fig. 1 and cropped at one FWHM. The two images correspond respectively to the (a) first and the (b) last frame.
In the text

Fig. 4. ROC curves for the NACO, SPHERE and LMIRCam data sets, with respectively, the KLIP-FM RSM map using the Gausssian maximum likelihood for the pre-optimisation of the flux parameter β (red), the KLIP-FM RSM map with no flux pre-optimisation (NF), which relies on the standard maximum likelihood used in the original RSM map for the estimation of flux parameter β (blue), the forward model matched filter KLIP-FMMF (green), the RSM map using KLIP (orange) and KLIP using the standard S/N map (black). (a) NACO at 2λ/D, (b) SPHERE at 2λ/D, (c) LMIRCam at 2λ/D, (d) NACO at 8λ/D, (e) SPHERE at 8λ/D, (f) LMIRCam at 8λ/D.

In the text

Fig. 5. ROC curves for the NACO, SPHERE and LMIRCam data sets, with, respectively, the LOCI-FM RSM using the Gausssian maximum likelihood for the pre-optimisation of the flux parameter β (red), the LOCI-FM RSM map with no flux pre-optimisation (NF), which relies on the standard maximum likelihood used in the original RSM map for the estimation of flux parameter β (blue), the RSM map using LOCI (orange) and LOCI using the standard S/N map (black). (a) NACO at 2λ/D, (b) SPHERE at 2λ/D, (c) LMIRCam at 2λ/D, (d) NACO at 8λ/D, (e) SPHERE at 8λ/D, (f) LMIRCam at 8λ/D.

In the text

	Fig. 6. ROC curves for the NACO, SPHERE and LMIRCam data sets, with the four best combinations of PSF subtraction techniques used to generate the RSM map algorithm. (a) NACO at 2λ/D, (b) SPHERE at 2λ/D, (c) LMIRCam at 2λ/D, (d) NACO at 8λ/D, (e) SPHERE at 8λ/D, (f) LMIRCam at 8λ/D.
In the text

	Fig. 7. Evolution of the probabilities for the forward, backward, and forward-backward approaches using KLIP, around the location of a planetary candidate injected in the 51 Eridani data set (radial distance of 4λ/D with a contrast of 3.76 × 10−5).
In the text

	Fig. 8. RMS position errors averaged over the two axes expressed in mas for the KLIP RSM map using the forward (blue) and forward-backward (red) versions of the RSM approach and for the NEGFC approach (orange). The two graphs show the dependence of the averaged RMS position error on the contrast for respectively a radial distance of (a) 2λ/D and (b) 8λ/D.
In the text

	Fig. 9. RSM map generated with forward and forward-backward approach for the LMIRCam data set using KLIP with 18 principal components, a FOV rotation expressed in terms of FWHM of 0.3 and a Gaussian distribution. A square root-based scale has been used to increase the background speckle noise brightness. (a) Original RSM map and (b) Forward-Backward RSM map.
In the text

Fig. 10. Procedure for the estimation of the contrast curve, the orange star corresponds to the linear interpolation, the black star to the corresponding TPR, the blue stars to the two previous estimations used for the linear interpolation and the green star to the flux corresponding to a TPR of 0.5. The red curve makes the link between the flux and the true positive rate in the case of a single false positive while the green line shows the tolerance’s interval of [0.45,0.55].

In the text

	Fig. 11. Contrast curves for the three data sets using the optimal combination of PSF subtraction techniques obtained in Sect. 5 (resp. 3, 11, and 10) with the region [2λ/D,16λ/D] considered for the SPHERE data set and the region [2λ/D,8λ/D] considered for the other two data sets to get the contrast in the first arcsecond. (a) NACO, (b) SPHERE, (c) LMIRCam.
In the text

Fig. B.1. ROC curves for the NACO, SPHERE, and LMIRCam data sets, with respectively the KLIP RSM map using the Gaussian maximum likelihood for the pre-optimisation of the flux parameter β (red), the KLIP RSM map with no flux pre-optimisation (NF), which relies on the standard maximum likelihood used in the original RSM map for the estimation of flux parameter β (blue), and the KLIP RSM map using the forward-backward approach for the probability estimation instead of the original forward approach (green). (a) NACO at 2λ/D, (b) SPHERE at 2λ/D, (c) LMIRCam at 2λ/D, (d) NACO at 8λ/D, (e) SPHERE at 8λ/D, (f) LMIRCam at 8λ/D.

In the text

	Fig. C.1. ROC curves for the NACO, SPHERE, and LMIRCam data sets, with the LOCI-FM RSM map using respectively a crop size for the froward modelled PSF of 5 (red), 7 (blue), 9 (green), 11 (orange) pixels (FWHM ≈ 5 pixels for all three data sets). (a) NACO at 2λ/D, (b) SPHERE at 2λ/D, (c) LMIRCam at 2λ/D, (d) NACO at 8λ/D, (e) SPHERE at 8λ/D, (f) LMIRCam at 8λ/D.
In the text

Fig. D.1. Cumulative maximum TPR reached without any FP (blue), and cumulative average number of FPs inside the entire frame at TPR = 1 (red) for 16 different combinations of PSF subtraction techniques used to generate the RSM map and two radial distances. The dark-coloured bars are obtained at a radial distance of 2λ/D while the light-coloured bars correspond to 8λ/D for, respectively, the NACO (top), SPHERE (middle), and LMIRCam (bottom) data sets. The asterisks for some values of FP at TPR = 1 indicate that a TPR of 1 has not been reached at a distance of 8λ/D and that the smallest probability threshold was chosen instead (highest TPR). A high performance for a combination of PSF subtraction techniques corresponds to a tall blue histogram alongside a short red histogram.

In the text

	Fig. E.1. Contrast curves for the SPHERE data set using KLIP RSM with 20 principal components and a FOV rotation expressed in terms of FWHM of 0.3. The region [2λ/D,16λ/D] has been considered to get the contrast in the first arcsecond.
In the text