Free Access
Issue
A&A
Volume 633, January 2020
Article Number A95
Number of page(s) 15
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201936421
Published online 17 January 2020

© ESO 2020

1 Introduction

High contrast imaging (HCI) is one of the most challenging techniques for exoplanet detection, but is also one of the most promising (see Bowler 2016, for a review). The main difficulties encountered with HCI arise from the small angular separation between the host star and the potential exoplanets, the flux ratio between them (usually below 10−3), and the image degradation caused by the Earth’s atmosphere. Adaptive optics (AO) and coronagraphic techniques are now widely used to improve the quality and reduce the dynamic range of the images in dedicated instruments such as GPI (Macintosh et al. 2008), SPHERE (Beuzit et al. 2019), or SCExAO (Lozi et al. 2018). However, despite the use of these cutting-edge technologies, the resulting images are still affected by residual aberrations. Under good observing conditions, the performance of HCI instruments is limited by aberrations arising in the optical train of the telescope and instrument, generating quasi-static speckles in the field of view. Different processing techniques along with observing strategies have been proposed in the last decade to deal with these quasi-static speckles, whose shape and intensity (about 10−4) are similar to potential companions.

Angular differential imaging (ADI, Marois et al. 2006) is nowadays the most commonly used observing strategy to mitigate quasi-static speckles in HCI. This observing strategy consists in acquiring images in pupil tracking mode, that is, with the instrument derotator keeping the pupil orientation fixed. The aim of this approach is to keep the quasi-static speckles fixed in the focal plane, so that they can easily be identified with respect to astrophysical objects rotating around the star along with the parallactic angle. Using this temporal diversity, a model of the speckle field, often referred to as the reference point spread function (PSF), may be built from the data. This reference PSF is then subtracted from the set of ADI images. The resulting residuals frames are eventually aligned and combined to detect the signal of potential exoplanets or discs, which should not have suffered too much from subtraction of the reference PSF.

Several methods using this approach have been proposed to maximise the noise reduction, such as for example a locally optimised combination of images (LOCI, Lafreniere et al. 2007), principal component analysis (PCA, Soummer et al. 2012), non-negative matrix factorization (NMF, Ren et al. 2018), and low rank plus sparse decomposition (LLSG, Gomez Gonzalez et al. 2016), allowing the user to reach contrasts down to 10−6 at 0.5 arcsec in the H-band (1.6 μm) with the latest generation of HCI instruments (e.g. Vigan et al. 2015). Another family of post-processing algorithms replaces the reference PSF subtraction by a forward modelling of the planetary companion using an inverse problem framework (ANDROMEDA, Cantalloube et al. 2015; FMMF, Pueyo 2016; Ruffio et al. 2017). In both cases, the detection is typically performed via the estimation of signal-to-noise-ratio (S/N) maps. In contrast with the forward- model-based algorithms, which provide a S/N map as a by-product of the model, in the case of reference PSF subtraction the S/N map is usually generated by using the median-averaged residual frames to estimate the annulus-wise S/N of every pixel it contains. The detection of planetary candidates is then done via the definition of an S/N threshold. Several methods have been proposed to generate S/N maps from the set of reference PSF-subtracted images (e.g. Mawet et al. 2014; Bottom et al. 2017; Pairet et al. 2019).

In this paper, we propose a novel approach to dealing with this last step of the ADI sequence post-processing. Instead of averaging the set of de-rotated images obtained after the reference PSF subtraction and computing an S/N map, we propose to consider the entire set of residual frames and rely on a regime-switching algorithm to classify the pixel values into two categories, regrouping either the planetary signals or the quasi-static speckles. The probability associated with the planetary regime then allows the creation of a detection map. The algorithm derives fromthe Markov regime-switching model first proposed by Hamilton (1988), which is widely applied to analyse economic and financial time series. The aim of our new detection algorithm is to more effectively treat the residual noise still observed in the cube of residuals provided by ADI methods, increasing our ability to disentangle faint signalsfrom bright speckles. The flexibility of the algorithm allows the use of ADI cubes treated with most post-processing methods. The cubes of residuals obtained from the different post-processing methods may be used separately, but can also be used together, further improving the sensitivity of the detection algorithm to faint companions.

The rest of the paper is organised as follows. In Sect. 2, we describe the new regime-switching model for the detection of exoplanets. Section 3 presents in detail the model estimation and the definition of the differentparameters. The ability of our model to disentangle faint planetary signals from bright speckles is tested in Sect. 4 by injecting fake companions into two different data sets and by comparing the results with state-of-the-art ADI-based post-processing techniques. Finally, Sect. 5 concludes on this work.

2 Regime-switching model

The proposed detection algorithm derives from the Markov-switching regressions introduced by Goldfeld & Quandt (1973) and Cosslett & Lee (1985) and further improved by Hamilton (1988, 1994), who developed an iterative inference algorithm to estimate themodel parameters, namely the Markov regime-switching model (RSM). This approach is one of the most popular non-linear time series models in the econometric literature and many variants have been proposed. The aim of the RSM is to take into account possible dramatic changes in the behaviour of time series such as the transition between economic expansion and contraction in the case of financial time series. The regime-switching model relies on several linear equations to describe the different states of a system described by a time series. The probability of being in a given state depends on both a pre-defined transition probability and on the ability of the different equations to properly describe the evolution of the time series. One of the model outcomes is the probability associated with different regimes. For each element of the time series, the RSM provides the probability of being in any one of the different regimes. Our detection map derives directly from these probabilities.

In the case of our RSM detection map, the time series is built from the de-rotated cube of residuals obtained after the PSF subtraction and de-rotation steps of the ADI sequence post-processing. Several cubes of residuals treated with different ADI PSF subtraction techniques may be stacked in the time axis to provide additional information and increase the ability of the model to detect faint companions. To allow for the detection of planetary signals, we rely on two different regimes to model the de-rotated cube of residuals: a regime in which the residuals time series is described by speckle noise and a second regime with speckle noise plus a planetary signal. The planetary signal may be modelled as the measured off-axis PSF1 or as a forward model of the off-axis PSF after the subtraction. We consider in this paper the measured off-axis PSF for simplicity, although the algorithm may be easily adapted to a forward modelled off-axis PSF.

The RSM we propose here is a modified version of the original Markov-switching model, in which only one parameter is determined via a maximum log-likelihood estimation. We rely on the characteristics of the data set to define the other model parameters. Having presented the basic principles behind our RSM, we may now describe the detailed procedure for our RSM detection map computation.

2.1 Building the time series

The first step of our estimation procedure is to build the time series that the regime switching model will try to model. As the noise properties are expected to evolve with radial distance, the regime switching model is applied annulus-wise. For each annulus a, a specific residuals time series is built by vectorizing that part of the cube of residuals, indexed by ia the flattened pixel number. The length of the time series depends on the number of pixels in the considered annulus La but also onthe number of frames in the original de-rotated cube of residuals T. We indeed take advantage of all the individual frames contained in the de-rotated cube of residuals instead of collapsing the cube as is usually done when estimating an S/N map. As can be seen from Fig. 1, the time series is built by concatenating the set of T observations for every pixel contained in the annulus a, i.e. with ia ∈ {1, …, T × La}. The first subscript of X indicates the selected frame in the de-rotated cube of residuals, while the second one provides the position of the considered pixel in the selected annulus a. Both subscripts are replaced by a single index ia to form the residuals time series that feeds the RSM.

We consider first the time axis and then the spatial axis in order to stay in the planetary regime during T steps of theiterative process used to build the detection map, instead of switching T times between both regimes when a planetary signal is present in a given annulus. Indeed, when travelling through the residuals time series, the planetary signal observed in a given pixel will act on the regime-switching model during T steps, allowing theprobability of being in the planetary regime to build up thanks to the short-term memory of the model. This helps to enhance the sensitivity of the algorithm to faint signals as it allows the probability to build up for a longer period of time.

2.2 Model description

The second step of the RSM detection map computation consists in defining the set of equations describing the residuals time series for the two considered regimes. In the first regime, the time series is described by a residual noise following the statistics of the quasi-static speckle residuals contained in the annulus. In the second regime, the time series is described by both the residual noise and the planetary signal model (off-axis PSF). The PSF being two-dimensional, we consider not only one pixel at a time but a batch of pixels in a square of size θ equal to the full width at half maximum (FWHM) of the PSF. In order to define the probability of observing a planetary signal at a given pixel , we therefore need to consider a number of neighbouring pixels depending on the value of θ. As depicted in Fig. 2, we define as the residuals matrices of dimension θ × θ centred on , which will replace the time series used so far. Larger values of θ may be considered in the case of a forward-modelled off-axis PSF to take into account the signal self-subtraction, which, for instance, could create negative wings in the azimuthal direction. Our RSM is therefore characterised by the following equations: (1)

where β provides the strength of the planetary signal, μ the mean of the quasi-static speckle residuals, and their time and space varying part characterised by the quasi-static speckle residuals statistics (see Table 1 for a summary of all the variables used in the RSM). Here, P is the model of the planetary signal, which is the normalised off-axis PSF in the FWHM region.

As can be seen from Eq. (1), there exist two possible states , which are reflected in the value taken by theparameter , with in the case of a planetary signal detection and in the other case. Here, is not directly observable, but we see its effect on the behaviour of via the realisation .

The parameter is a realisation of a two-state Markov chain allowing short-term memory. This implies that we only consider the state in which the system was at index ia − 1 to define the probability of being in agiven state for the current index ia. The fact that therealisation is a probabilistic outcome implies that we cannot consider being in only one of the two regimes. We have instead a given probability of being in each of them. Our RSM tries to describe the behaviour of the time series via a probability-weighted sum of the values generated by the equation describing each regime.

thumbnail Fig. 1

Residuals time series for a given annulus a is obtained by stacking the pixel values of the considered annulus along the time axis.

Open with DEXTER
thumbnail Fig. 2

Residuals matrices obtained from the first frame of the cube of residuals for the last three pixels of the annulus with θ equal to 3. The time series is created by considering matrices of dimension θ × θ centred on every in the cube of residuals.

Open with DEXTER

2.3 Definition of the model probabilities

The probability of being in a state or regime is characterised by the set of parameters of Eq. (1), that is, P the planetary signal model, and μ and β, the statistical properties of the residual noise . We make the simplifying assumption here that the quasi-static speckles residuals may be characterised to a good level of precision by their mean μ and variance σ. We write the probability of observing in the state s at step ia as follows: (2)

where P, μ, β, σ and provide the parameters of the model.

This probability is the key element of our RSM detection map as the map is constructed based on the value taken by for every pixel of every annulus. Indeed, provides a detection probability for each pixel and each frame of the de-rotated cube of residuals. The final RSM detection map is created by averaging these probabilities along the time axis of the cube of residuals.

In the case of a two-state Markov chain, the computation of necessitates the estimation of (i) the probability of observing the system in the state q at step ia − 1, (ii) the transition probability pq,s from state q to state p and (iii) the likelihood of observing in state s at step ia, which we note . The probability of being in a state s at index ia can be computed as the normalised likelihood of being in state s at index ia multiplied by the probability of having been in either of the two states at index ia −1 and by the transition probability pq,s, which accounts for the short-term memory of the algorithm. The expression of the state probability is therefore given by the following expression (Hamilton 1988): (3)

with the sum f of conditional densities for index ia given by: (4)

and the transition probabilities given by: (5)

with q, s ∈{0, 1}. We consider the two possible states describing the system at index ia − 1 via the sum over q. The function , which represents the numerator summed over the two possible states taken at index ia, ensures that the sum of the probability equals one for every index ia.

Table 1

Description of the mathematical notations for the variables used in the RSM detection map computation.

2.4 Transition probabilities estimation

For our two-regime model, the transition probability pq,s regroups the probabilities of staying in either regime along with the probabilities of switching to the other regime. The estimation of pq,s is relatively straightforward by imposing on the algorithm the potential existence of no more than one planetary signal per annulus. A number of planetary signals per annulus in the interval ]0,1] may therefore be considered. Following our testing, a value of one companion per annulus must be privileged in the case of faint companions as lower values decrease both the residual speckles and the companion intensities in our model. Considering the number of pixels La and the number of frames T, the parametrisation of pq,s translates as follows in the case of one planetary signal per annulus: (6)

2.5 Likelihood function definition

The determination of the likelihood is the key step of the model estimation. The challenge is to select the right probability distribution function to properly describe , the residual noise due to the quasi-static speckles. Indeed, the value taken by depends directly on the position of the elements of , or the elements of , in the probability distribution of the quasi-static speckle residuals. Considering the small transition probabilities p0,1, the probability of planetary signal detection depends heavily on the value taken by . The parametrisation of the selected probability distribution function also plays an important role.

Different probability distribution functions may be used. For the sake of clarity, we illustrate the likelihood function definition with a simple Gaussian distribution as is done in Hamilton (1988). However, the following section will allow us to investigate the question of the optimal probability distribution function selection as different post-processing algorithms provide different noise distributions for different separations. The Gaussian distribution allows us to construct a likelihood function for state s at index ia in the following manner: (7)

with n the index of the matrix elements for and P. The sum over the matrix elements allows us to obtain only one value per considered θ × θ patch.

2.6 Model estimation

Since the estimation of depends on its value at the previous step, we rely on an iterative procedure to estimate the entire set of . This iterative procedure requires the definition of an initial condition for ξq,0. Assuming that the considered Markov chain is ergodic, we can simply set ξq,0 = P(St = qP, μ, β, σ) equal to the unconditional probability ξq,0 = P(St = q). Following the approach proposed by Hamilton (1994), the two initial probabilities ξ0,0 and ξ1,0 may be estimated using the following system of equations: (8)

which translates in terms of matrices into: (9)

with the set of initial probabilities, , and A given by: (10)

with P the matrix of pq,s, I2×2 a diagonal matrix of dimension 2 × 2. Solving the system of Eq. (8) to obtain the initial probabilities, ξ, is then equivalent to taking the third row of the matrix .

3 Detection map estimation

We propose in this section a procedure to produce a RSM detection map. The model we developed so far necessitates the computation of cubes of residuals along with the definition of several parameters: the probability distribution function of the quasi-static speckles residuals and its first two moments, the planetary signal model P, the intensity parameter β, and the transition probability pq,r. The transition probability pq,r is already defined in Sect. 2.4. We therefore consider the remaining three model parameters.

3.1 Computation of de-rotated cubes of residuals

The first step to create a RSM detection map is the production of the de-rotated cubes of residuals for the selected ADI-based post-processing techniques feeding our regime-switching algorithm. As an illustration of the ability of our model to improve the detection when considering several methods at once, in this paper we consider three different post-processing techniques: annular PCA, NMF, and LLSG. For the two first approaches, the estimation of the cubes of residuals starts with the definition of a reference PSF. Annular PCA follows the PCA principles by computing the directions of maximal variance from the main matrix representing the ADI sequence, , with n the number of frames and p the number of pixels in the considered annulus. The determination of a reference PSF is done via the estimation of the eigenvectors V of the matrix M by taking Vk, the first k components of V. Annular PCA relies on a separate estimation for each annulus composing the original cube of data to take into account the radial evolution of the noise distribution; it allows the user to consider the local structure of the speckle noise instead of the entire frame. The cube of residuals is then obtained via the subtraction of the low rank matrix from the initial ADI sequence M.

As for annular PCA, NMF can be understood as a low rank approximation, with an additional non-negativity condition. This method consists in the decomposition of a matrix into two factors of non-negative values via the minimisation of the Frobenius norm: (11)

where and . The method allows the definition of a matrix WH with rank k lower than that of the original matrix M, keeping only the main components of M. The matrix WH provides a reference PSF for the entire set of frames representing the structure of the residual starlight. As for annular PCA, this matrix is subtracted from the original ADI sequence to obtain the cube of residuals, MWH.

Finally, the LLSG estimation is based on the decomposition of the cube intensities in three separate components: L, a low-rank matrix, S, a sparse matrix expected to contain the potential planetary signal, and G, the Gaussian part of the background noise. This partly explains why the distribution of the resulting residuals observed in Fig. 3 is far from being Gaussian, the Gaussian part of the noise having already been removed. More information about the algorithm may be found in Gomez Gonzalez et al. (2016). The cube of residuals is directly provided by S.

3.2 Probability distribution function

We then move to the model parameters definition by first considering the selection of the probability distribution function describing the speckle residuals. Figures 3a–d provides the distribution of the residuals for a VLT/NACO ADI sequence (see Sect. 4 for a description of the data set) obtained with respectively the annular PCA, the NMF, and the LLSG methods. We see from these graphs that the distribution of the residuals is either close to a Lapacian or to a Gaussian distribution depending on the selected post-processing techniques and on the angular separation. At small angular separations, the tails of the distributions of the residuals seem to be closer to a Laplacian, while at larger separation they seem closer to a Gaussian, except for LLSG processing. This radial evolution is mainly due to the higher (relative) number of intense speckles, the lower number of pixels, and the lower field rotation at small separation. Overall, the distribution of the residuals is close to a Gaussian for annular PCA and NMF, and close to a Laplacian for LLSG. This partially confirms the findings of Pairet et al. (2019), who demonstrated that the residuals were closer to a Laplacian than a Gaussian distribution, especially when looking at the tails of the distribution.

The results of Fig. 3 illustrate the difficulty of defining the residuals distribution as there exists a dependence on both the separation and the post-processing technique along with differences between the tails and the core of the distribution. We therefore consider both the Gaussian and Laplacian distributions in the performance assessment of Sect. 4.

The proposed regime-switching model provides a local detection probability as it considers one annulus at a time. The parameters of the residuals probability distribution should therefore be estimated locally. In the previous section, we considered not a single pixel at a time but a θ × θ matrix of pixels centred on the pixel of interest. We therefore estimate the pixel-wise mean and variance of the residuals empirically by considering an annulus with a width of θ pixels centred on the selected annulus. The entire set of frames is used for the estimation of these two parameters. Although planetary signal may be included in the annulus, the effect of this signal on the estimation of the mean and variance is limited and decreases with angular separation.

3.3 Intensity parameter

For the estimation of the intensity parameter β, we rely on the estimated variance of the pixel intensity in the annulus. We were inspired here by the S/N maps that are usually created with the final frame provided by most of the ADI techniques. We define the intensity parameter β as a multiple of the estimated variance σ: (12)

The β parameter is the only parameter we propose to estimate via a maximum log likelihood. Several values of δ are tested in a given interval starting at δ = 1, as δ = 0 would imply a single regime model. The optimal δ in an annulus a is the one leading to the highest log-likelihood sum in the considered interval.

Relying on this definition of β allows us to obtain information about the position of the detected planetary signal inside the probability distribution of the residual speckles. A higher δ implies that the detected signal is farther in the distribution tails, which indicates a higher level of confidence (which will generally translate into a higher probability in the RSM map) about the detected planet and a higher flux for a given noise distribution. However, the β parameter does not provide an estimation of the planetary flux as we are not using a forward model of the PSF for the planetary signal.

thumbnail Fig. 3

Distribution of the residuals for a VLT/NACO data set after PSF subtraction by annular PCA (top), NMF (middle), and LLSG (bottom) along with a Gaussian (orange line) and Laplacian (green line) fit at small (left) and large separations (right), with respectively 20 components for the annular PCA and the NMF and a rank of 5 for the LLSG.

Open with DEXTER

3.4 Planetary signal model

Using a forward model for the planetary signal would allow us to take into account the distortions (such as self-subtraction) created by ADI-based post-processing treatment when estimating cubes of residuals. Although a forward-modelled PSF should provide more accurate results, it should be noted that some ADI PSF subtraction techniques do not lend themselves to the analytical computation of a forward model (e.g. LLSG, NMF). A more universal numerical way to compute a forward model is to compare theinitial cube of residuals and the one in which a fake companion has been injected. Following this approach, we tested numerical estimation of forward modelled PSF for generating RSM detection maps but without managing to improve the algorithm accuracy compared to the use of measured off-axis PSF. We therefore decided to only consider measured off-axis PSFs in the rest of this paper. However, a forward model variant of the proposed algorithm is still under development and should be a valuable improvement of the current model, at the expense of its computation time.

3.5 Regime switching model detection map estimation

Now that we have defined the procedure to estimate the cubes of residuals feeding the RSM algorithm as well as the model parameters, we may summarise the main steps of the algorithm as follows:

  • 1.

    Compute the residuals cubes for the selected ADI techniques and de-rotate all the resulting frames;

  • 2.

    define the separation to the star for the first and last annuli, respectively aini = FWHM∕2 + 1 and afin = (fsizeFWHM)∕2 with fsize the size of the frame;

  • 3.

    define the series for the first annulus;

  • 4.

    estimate the mean and variance of the residuals inside the annulus separately for each residuals cube;

  • 5.

    using the iterative procedure described in Sect. 2, estimate for each index ia for the set of tested δ;

  • 6.

    include the probability of planetary signal providing the maximum likelihood in a three-dimensional matrix ;

  • 7.

    repeat steps 2–6 for the next annulus (a + 1) until afin is reached;

  • 8.

    average the detection probability contained in U along the time axis to obtain the final RSM detection map.

The resulting detection map provides the averaged probability of observing a planetary signal in a given cube of data, along with the optimal β. The following section explores the effectiveness of this new approach when applied to observational data sets.

4 Performance assessment

4.1 Data

We propose the use of two ADI sequences acquired with two instruments of the Very Large Telescope (VLT): NACO and SPHERE. This allows us to investigate the ability of our model to deal with the different noise profiles produced by these instruments.

The first data set focuses on β Pictoris and its planetary companion β Pictoris b. Itwas 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 averaging 40 successive individual exposures, each frame providing an effective integration time of 8 s. The parallactic angle ranges from −15° to +68°. We use every third frame to reduce the CPU time and cropped the central 101 × 101 pixels region to consider mainly the first arc-second.

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

4.2 Detection maps

We start our analysis by considering the RSM detection map generated with the proposed algorithm and based on the residual cubes provided by annular PCA, NMF, and LLSG, and compare it with the S/N map obtained with the same three post-processing algorithms. The post-processing as well as the S/N detection maps are generated for all three methods with the VIP package developed by Gomez Gonzalez et al. (2017) using the standard parametrisation. Both annular PCA and LLSG are performed annulus-wise, with each annulus being divided into four segments in the case of LLSG. Other parametrisations are possible as the proposed approach works with any de-rotated cube of residuals. The three cubes of residuals obtained with the selected post-processing techniques are then stacked to create a single cube to feed the RSM. The variance and the mean of the residuals are estimated separately for each subcube as their noise profiles are specific, as demonstrated in the previous section when looking at the residual distributions.

Figure 4 displays the RSM detection map and the S/N maps obtained for the SPHERE-IRDIS 51 Eridani data set (see Fig. A.1 for similar detection maps for the NACO β Pictoris data set). As can be seen, the difference in intensity between the planetary signal and the background speckles is much higher with our new approach than with the usual S/N maps. Fifty-one Eridani b (contrast of 6.73 × 10−6 ± 9.02 × 10−7 at a separation of 453.4 ± 4.6 mas, Samland et al. 2017; Maire et al. 2019) can be clearly identified on the lower left quadrant with RSM, annular PCA, and LLSG, although we observe a higher number of false positives in the case of LLSG. The visual identification becomes more difficult when looking at the S/N map provided by NMF, which shows brighter wind-driven halo residuals.

To illustrate the computation of the RSM map, Fig. 5 shows how the probability builds up when getting closer to a planetary signal; it reports the RSM map probabilities along the radial axis crossing the peak value attributed to 51 Eridani b, along with the optimal δ for the respective annuli. The data includes 7 pixels × 197 frames × 3 ADI-based post-processing techniques and is centred on the pixel showing the highest probability. As can be seen, no signal may be found in the first 591 patches representing the first pixel. The probability builds up steadily for the next three pixels until reaching a peak probability of over 95%. The value of the optimal δ increases as well with a peak value of 4 reached at the fifth pixel, illustrating the displacement of the signal farther into the residuals distribution tail due to the increasing flux coming from the planetary candidate. We then observe a decrease of the probability and optimal δ, which eventually gets back to the background speckle noise level. The stacked cube of residuals encompasses the cubes of residuals generated first by the annular PCA, then by the NMF, and finally by the LLSG. Looking at the sharp increase observed at the beginning of every pixel, we seethat the strongest signal may be found in the annular PCA cube of residuals, confirming the visual analysis of the S/N maps. However thesignal is still strong in the two other cubes of residuals to be able to maintain the high probability observed for the three central pixels. Changing the order of the cubes of residuals when computing the RSM detection map does not significantly affect the probabilities.

thumbnail Fig. 4

Probability map obtained for the SPHERE-IRDIS 51 Eridani data set, with the RSM using a Gaussian distribution along the S/N map generated with the cube of residuals obtained with annular PCA, LLSG, and NMF. The annular PCA and the NMF use 20 components, and the LLSG has a rank of 7. The colour scale indicates the probability for the RSM map and the S/N for the three S/N maps (Mawet et al. 2014). The maps are centred on the star 51 Eridani while 51 Eridani b is identified by the white circle in the lower left quadrant.

Open with DEXTER

4.3 Receiver operating characteristic curves

In order to explore in more detail the properties of our new approach and compare its performance with other state-of-the-art methods, we generated synthetic data sets based on the two ADI sequences presented in the previous section. We rely on the injection of fake companions in the initial ADI sequences, an approach widely accepted by the HCI community for generating synthetic data to assess the sensitivity of post-processing methods. Since the contrast that can be reached as well as the noise structure both depend on the angular separation, we consider three different annuli as described in Table 2. The comparison with the other methods is based on receiver operating characteristic (ROC) curves, which are widely used to assess the performance of binary classifiers. In these curves, one axis provides the true positive rate and the other the false positive rate. When using ROC curves for performance assessment, the main proxy for the classifier performance is the area under the ROC curve: the better the classifier, the higher the area under the ROC curve, i.e. the higher the true positive rate for a given false positive rate. We replace the false positive rate by the number of false positives for the entire frame, averaged over the number of test data sets considered for a given separation as is done in Gomez Gonzalez et al. (2018).

The fake companions are defined as the normalised off-axis PSF, generally measured by offsetting the target star from the coronograph, multiplied by flux values from a predefined interval defined to challenge the set of tested methods. Five different flux values are tested for each separation with step size of 0.5 times the initial value. For each flux value, eight positions are tested to mitigate the impact of bright speckles or local minima. The resulting 40 test data sets are then used to estimate the ROC curves for each separation. The contrasts for the three selected separations are provided for the NACO and SPHERE data sets in Table 2. Before injecting the fake companion, we removed the known companions and some bright disc structures for the β Pictoris data set using the negative fake companion technique (Lagrange et al. 2010). We consider a false positive to be a detected companion at any other location than the one selected for the fake companion injection.

The exoplanet detections for the annular PCA, the NMF, and the LLSG methods are based on S/N maps generated using the procedure of Mawet et al. (2014).The detection of a true or false positive is done on the de-rotated median-combined individual frame by estimating the S/N for every pixel. This estimation is done annulus-wise in order to take into account the evolution of the residuals distribution. The S/N is calculated by comparing the flux inside an aperture with a diameter of one FWHM centred on the considered pixel (i.e., 5 pixels for both data sets) with the flux of all the other apertures included in the annulus (for more details about the estimation see Mawet et al. 2014). Once the S/N map is computed, successive thresholds are applied onto the S/N map to create the ROC curves. For each threshold, the detection of the fake companion as well as the number of false positives are recorded and averaged over the entire set of synthetic data sets generated for the considered annulus to construct our false and true positive rates. We follow a similar procedure for the RSM detection map, simply replacing the S/N thresholds by successive percentage thresholds applied to the detection map.

The parameters of the different post-processing techniques have been selected to maximise the area under the ROC curves, that is, to maximise the true positive rate while minimising the number of false positives. For annular PCA and NMF, the number of principal components used to construct the reference PSF was set to 20 for both data sets. As for LLSG, we selected a rank value of 5 for the estimation of the matrix S for the β NACO data set and 7 for the SPHERE-IRDIS data set. As regards the RSM, the mean and variance of the residuals distribution are again estimated annulus-wise. As the fake companions injected into our simulations have relatively low flux values, we tested δ in the interval [1,5] and kept the one leading to the highest total log-likelihood to generate the final RSM map.

As an illustration of the detection map calculation for the generation of ROC curves, Fig. 6 shows the probability and S/N maps obtained by injecting fake companions with high contrast values at three different separations from the star 51 Eridani (2, 4 and 8 λD). As can be seen, apart from the signal injected at 8 λD which appears relatively clearly in the S/N for all three post-processing methods, the RSM map is the only map providing a clear detection for all three fake companions. A set of detection maps is shown in Fig. A.2 for the NACO β Pictoris data set, leading to similar conclusions.

thumbnail Fig. 5

Evolution of the probability in the RSM detection map around the location where 51 Eridani is detected along with the optimal δ for the respective annuli.

Open with DEXTER
Table 2

Injected companions contrast range for the three considered separations.

4.3.1 Influence of the probability distribution

We now turn to the estimation of the ROC curves which will provide more comprehensive results. We start by considering two different variants of RSM to investigate the choice of the probability distribution for the likelihood function definition. The two variants presented in Fig. 7 use the Gaussian and Laplacian distribution, respectively, to construct the likelihood function appearing in . The ROC curves are estimated for different separations; as we have seen in the previous section, the probability distribution describing the residuals evolves withangular separation. As can be seen from Fig. 7, the results of the two variants are very close in the case of the β Pictoris data set, while the distance between them becomes significant for the 51 Eridani data set. In both cases, the RSM model using the Laplacian distribution performs better for small separation while the Gaussian distribution leads to better results for larger separations.

These results confirm the findings made with Fig. 3 and the importance of tails fit when selecting the optimal probability distribution. It demonstrates the interest of considering the residuals distribution evolution along the radial axis to optimally parametrise our model. We therefore propose to start the RSM detection map estimation with an analysis of the noise profile to select the right probability distribution for every separation. This additional step has been included in the RSM detection map python package that we developed based on the model presented in this paper2. The function allows the user (i) to select one of the two distributions, (ii) to automatically select the best distribution based on a best-fit approach, or (iii) to create a hybrid distribution consisting in a weighted sum of both distributions. This last possibility can be useful when facing asymmetrical probability distributions as the parameters of both distributions may be estimated separately based on a best-fit approach.

4.3.2 Comparison with S/N-based detection

We now address the question of the performance of our algorithm compared to the three post-processing methods using S/N maps. For the two data sets, Fig. 8 reports the ROC curves of all four methods for the same separations as before. Considering the results presented in Fig. 7, we selected for each data set and each separation the distribution that provided the highest area under the ROC curve. The results demonstrate the interest of the new approach considering that the RSM performs better in every case. This may be explained by the ability of our model to be fed with multiple cubes of residuals, but also by its ability to focus only on relevant data thanks to the regime-switching feature. This allows our model to take advantage of the strength of the different post-processing methods used to produce the cubes of residuals. As speckles are not treated equally by these post-processing techniques, it is easier to remove them by taking into account several cubes of residuals. This ability to remove speckles is further improved by the memory of the RSM. Indeed, the dependence of on the transition matrix pq,s and on the probabilitiesat step ia − 1 (see Eq. (3)) partly mitigates the effect of speckles on the detection map. Outliers caused by quasi-static speckles do not lead to a clear regime switch, while when facing a planetary signal the detection probability builds up along the time axis as we see in Fig. 5. The dependence on the past observation reduces the noise in the final detection map significantly.

Furthermore, the possibility of selecting the right probability distribution to describe the residuals allows us to more precisely describe the behaviour of these residual speckles, which is not possible with the S/N approach. The more significant improvements for the 51 Eridani data set may be explained by the lower level of noise inside this ADI sequence, which suggests that our model should perform better with the latest generation of instruments.

thumbnail Fig. 6

Detection map obtained after injecting three fake companions in the SPHERE-IRDIS 51 Eridani reference cube used for the ROC estimation, at a distance of 2, 4, and 8 λD with a contrast of 1.0 × 10−4, 1.2 × 10−5 and 3.7 × 10−6, respectively.The colour scale indicates the probability for the RSM map and the S/N for the three S/N maps. The maps are centred on the star 51 Eridani while the fake companions are identified by the white circles.

Open with DEXTER

5 Conclusion

Here, we explore the possibility of improving exoplanet detection using an RSM derived from the field of econometrics, with one regime representing the planetary signal in addition to the speckle noise and the other only the speckle noise. This novel approach allows the creation of probability maps based on cubes of residuals obtained with different ADI-based post-processing techniques. The RSM algorithm can be associated with any ADI-based post-processing techniques as it can be fed with different cubes of residuals separately or jointly. The short memory process at the heart of our RSM detection map allows quasi-static speckles to be treated more effectively when using several cubes of residuals provided by different post-processing algorithms and thereby allows the user to reach better detection performance.

The RSM is easy to use as most of the parameters are estimated empirically. The only parameter that may need to be tuned is δ, which defines the strengthof the signal coming from the planetary candidates. The model automatically selects this parameter via a maximum log-likelihood approach. However, an upper value has to be defined for the interval. The estimation of the RSM map takes between three and ten times longer than the standard Mawet et al. (2014) S/N map computation time, depending on the size of the ADI sequence and on the upper value for the parameter δ.

We demonstrate the interest of our approach by injecting fake companions into two data sets provided by the VLT/NACO and VLT/SPHERE instruments. We compared the proposed RSM map with standard S/N maps obtained with three state-of-the-art methods: annular PCA, NMF, and LLSG. The ROC curves clearly demonstrate the interest of our model as it outperformsall the other methods for the three angular separations considered, and for both data sets. The results also confirm that the probability distribution of the residuals evolves with radial distance and that it should be taken into account in our model when defining the likelihood function used to estimate the probability of being in one of the two regimes. Indeed, the Laplacian distribution clearly performs better for close separations while the Gaussian one provides better results for larger angular distances. The possibility of optimally selecting the probability distribution based on the residual noise profile has been included in the RSM detection map python package that we have developed.

thumbnail Fig. 7

Receiver operating characteristic curves for the β Pictoris and 51 Eridani data sets, with the RSM using a Gaussian (blue) and Laplacian (red) distribution, respectively, to construct the likelihood function.

Open with DEXTER
thumbnail Fig. 8

ROC curves for the β Pictoris and 51 Eridani data sets. The RSM variant (in blue) providing the highest area under the ROC curve in Fig.7 has been selected. The red, yellow, and green ROC curves are computed using the S/N map generated with respectively the annular PCA, the NMF and the LLSG.

Open with DEXTER

Acknowledgements

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 colleague A.-L. Maire for sharing the SPHERE data set.

Appendix A NACO β Pictoris

thumbnail Fig. A.1

Probability map obtained for the NACO β Pictoris data set, with the RSM using a Gaussian distribution along the S/N map generated with the cube of residuals obtained with annular PCA, LLSG, and NMF. The annular PCA and the NMF use 20 components and the LLSG has a rank of 5. The colour scale indicates the probability for the RSM map and the S/N for the three S/N maps. The maps are centred on the star β Pictoris while β Pictoris b is identified by the white circle in the lower left quadrant.

Open with DEXTER
thumbnail Fig. A.2

Detection map obtained after injecting three fake companions in the NACO β Pictoris reference cube used for the ROC estimation at a distance of 2, 4, and 8 λD with a contrast of 3.3 × 10−4, 0.4 × 10−4 and 1.7 × 10−5, respectively. The colour scale indicates the probability for the RSM map and the S/N for the three S/N maps. The maps are centred on the star β Pictoris while the fake companions are identified by the white circles.

Open with DEXTER

References

  1. Absil, O., Milli, J., Mawet, D., et al. 2013, A&A, 559, L12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  2. Beuzit, J.-L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  3. Bottom, M., Ruane, G., & Mawet, D. 2017, Res. Notes AAS, 1, 30 [NASA ADS] [CrossRef] [Google Scholar]
  4. Bowler, B. P. 2016, PASP, 128, 102001 [NASA ADS] [CrossRef] [Google Scholar]
  5. Cantalloube, F., Mouillet, D., Mugnier, L. M., et al. 2015, A&A, 582, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Cosslett, S. R., & Lee, L.-F. 1985, J. Econom., 27, 79 [Google Scholar]
  7. Delorme, P., Meunier, N., Albert, D., et al. 2017, SF2A-2017: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, 347 [Google Scholar]
  8. Goldfeld, S. M., & Quandt, R. E. 1973, J. Econom., 1, 3 [Google Scholar]
  9. Gomez Gonzalez, C., Absil, O., Absil, P.-A., et al. 2016, A&A, 589, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Gomez Gonzalez, C. A., Wertz, O., Absil, O., et al. 2017, AJ, 154, 7 [NASA ADS] [CrossRef] [Google Scholar]
  11. Gomez Gonzalez, C. A., Absil, O., & Van Droogenbroeck M. 2018, A&A, 613, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Hamilton, J. D. 1988, J. Econ. Dyn. Control, 12, 385 [Google Scholar]
  13. Hamilton, J. D. 1994, Time Series Analysis (Princeton: Princeton University Press) [Google Scholar]
  14. Lafreniere, D., Marois, C., Doyon, R., Nadeau, D., & Artigau, E. 2007, ApJ, 660, 1 [Google Scholar]
  15. Lagrange, A.-M., Bonnefoy, M., Chauvin, G., et al. 2010, Science, 329, 57 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  16. Lozi, J., Guyon, O., Jovanovic, N., et al. 2018, Adaptive Optics Systems VI (Bellingham, USA: SPIE Press) [Google Scholar]
  17. Macintosh, B. A., Graham, J. R., & Palmer, D. W. 2008, SPIE Conf. Ser., 7015, 701518 [Google Scholar]
  18. Maire, A.-L., Rodet, L., Cantalloube, F., et al. 2019, A&A, 624, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Marois, C., Lafreniere, D., Doyon, R., Macintosh, B., & Nadeau, D. 2006, ApJ, 641, 1 [Google Scholar]
  20. Mawet, D., Milli, J., Wahhaj, Z., et al. 2014, ApJ, 792, 97 [NASA ADS] [CrossRef] [Google Scholar]
  21. Pairet, B., Cantalloube, F., Gomez Gonzalez, C., Absil, O., & Jacques, L. 2019, MNRAS, 487, 2262 [NASA ADS] [CrossRef] [Google Scholar]
  22. Pueyo, L. 2016, ApJ, 824, 117 [NASA ADS] [CrossRef] [Google Scholar]
  23. Ren, B., Pueyo, L., Zhu, G. B., Debes, J., & Duchêne, G. 2018, ApJ, 852 [Google Scholar]
  24. Ruffio, J.-B., Macintosh, B., Wang, J. J., & Pueyo, L. 2017, ApJ, 842, 1 [Google Scholar]
  25. Samland, M., Mollière, P., Bonnefoy, M., et al. 2017, A&A, 603, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Soummer, R., Pueyo, L., & Larkin, J. 2012, ApJ, 755, L2 [Google Scholar]
  27. Vigan, A., Gry, C., Salter, G., et al. 2015, MNRAS, 454, 129 [NASA ADS] [CrossRef] [Google Scholar]

1

For coronagraphic imaging, an off-axis non-coronagraphic image of the target is routinely acquired before and after the observing sequence. This PSF reference is used to calibrate the flux of the star and provide a model of the planetary signal for a forward model-based algorithm. For non-coronagraphic imaging, this reference PSF is the unsaturated exposure.

2

The RSM detection map python package is available on GitHub: https://github.com/chdahlqvist/RSMmap

All Tables

Table 1

Description of the mathematical notations for the variables used in the RSM detection map computation.

Table 2

Injected companions contrast range for the three considered separations.

All Figures

thumbnail Fig. 1

Residuals time series for a given annulus a is obtained by stacking the pixel values of the considered annulus along the time axis.

Open with DEXTER
In the text
thumbnail Fig. 2

Residuals matrices obtained from the first frame of the cube of residuals for the last three pixels of the annulus with θ equal to 3. The time series is created by considering matrices of dimension θ × θ centred on every in the cube of residuals.

Open with DEXTER
In the text
thumbnail Fig. 3

Distribution of the residuals for a VLT/NACO data set after PSF subtraction by annular PCA (top), NMF (middle), and LLSG (bottom) along with a Gaussian (orange line) and Laplacian (green line) fit at small (left) and large separations (right), with respectively 20 components for the annular PCA and the NMF and a rank of 5 for the LLSG.

Open with DEXTER
In the text
thumbnail Fig. 4

Probability map obtained for the SPHERE-IRDIS 51 Eridani data set, with the RSM using a Gaussian distribution along the S/N map generated with the cube of residuals obtained with annular PCA, LLSG, and NMF. The annular PCA and the NMF use 20 components, and the LLSG has a rank of 7. The colour scale indicates the probability for the RSM map and the S/N for the three S/N maps (Mawet et al. 2014). The maps are centred on the star 51 Eridani while 51 Eridani b is identified by the white circle in the lower left quadrant.

Open with DEXTER
In the text
thumbnail Fig. 5

Evolution of the probability in the RSM detection map around the location where 51 Eridani is detected along with the optimal δ for the respective annuli.

Open with DEXTER
In the text
thumbnail Fig. 6

Detection map obtained after injecting three fake companions in the SPHERE-IRDIS 51 Eridani reference cube used for the ROC estimation, at a distance of 2, 4, and 8 λD with a contrast of 1.0 × 10−4, 1.2 × 10−5 and 3.7 × 10−6, respectively.The colour scale indicates the probability for the RSM map and the S/N for the three S/N maps. The maps are centred on the star 51 Eridani while the fake companions are identified by the white circles.

Open with DEXTER
In the text
thumbnail Fig. 7

Receiver operating characteristic curves for the β Pictoris and 51 Eridani data sets, with the RSM using a Gaussian (blue) and Laplacian (red) distribution, respectively, to construct the likelihood function.

Open with DEXTER
In the text
thumbnail Fig. 8

ROC curves for the β Pictoris and 51 Eridani data sets. The RSM variant (in blue) providing the highest area under the ROC curve in Fig.7 has been selected. The red, yellow, and green ROC curves are computed using the S/N map generated with respectively the annular PCA, the NMF and the LLSG.

Open with DEXTER
In the text
thumbnail Fig. A.1

Probability map obtained for the NACO β Pictoris data set, with the RSM using a Gaussian distribution along the S/N map generated with the cube of residuals obtained with annular PCA, LLSG, and NMF. The annular PCA and the NMF use 20 components and the LLSG has a rank of 5. The colour scale indicates the probability for the RSM map and the S/N for the three S/N maps. The maps are centred on the star β Pictoris while β Pictoris b is identified by the white circle in the lower left quadrant.

Open with DEXTER
In the text
thumbnail Fig. A.2

Detection map obtained after injecting three fake companions in the NACO β Pictoris reference cube used for the ROC estimation at a distance of 2, 4, and 8 λD with a contrast of 3.3 × 10−4, 0.4 × 10−4 and 1.7 × 10−5, respectively. The colour scale indicates the probability for the RSM map and the S/N for the three S/N maps. The maps are centred on the star β Pictoris while the fake companions are identified by the white circles.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.