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/00046361/201936421  
Published online  17 January 2020 
Regimeswitching model detection map for direct exoplanet detection in ADI sequences
^{1}
STAR Institute, Université de Liège,
Allée du Six Août 19c,
4000
Liège,
Belgium
email: carlhenrik.dahlqvist@uliege.be
^{2}
Max Planck Institute for Astronomy,
Königstuhl 17,
69117
Heidelberg,
Germany
Received:
31
July
2019
Accepted:
5
December
2019
Context. Beyond the choice of wavefront control systems or coronographs, advanced data processing methods play a crucial role in disentangling potential planetary signals from bright quasistatic speckles. Among these methods, angular differential imaging (ADI) for data sets obtained in pupil tracking mode (ADI sequences) is one of the foremost research avenues, considering the many observing programs performed with ADIbased techniques and the associated discoveries.
Aims. Inspired by the field of econometrics, here we propose a new detection algorithm for ADI sequences, deriving from the regimeswitching model first proposed in the 1980s.
Methods. The proposed model is very versatile as it allows the use of PSFsubtracted data sets (residual cubes) provided by various ADIbased techniques, separately or together, to provide a single detection map. The temporal structure of the residual cubes is used for the detection as the model is fed with a concatenated series of pixelwise time sequences. The algorithm provides a detection probability map by considering two possible regimes for concentric annuli, the first one accounting for the residual noise and the second one for the planetary signal in addition to the residual noise.
Results. The algorithm performance is tested on data sets from two instruments, VLT/NACO and VLT/SPHERE. The results show an overall better performance in the receiver operating characteristic space when compared with standard signaltonoiseratio maps for several stateoftheart ADIbased postprocessing algorithms.
Key words: techniques: image processing / techniques: high angular resolution / methods: statistical / methods: data analysis / planetary systems / planets and satellites: detection
© 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 cuttingedge 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 quasistatic speckles in the field of view. Different processing techniques along with observing strategies have been proposed in the last decade to deal with these quasistatic 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 quasistatic 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 quasistatic 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), nonnegative 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 Hband (1.6 μm) with the latest generation of HCI instruments (e.g. Vigan et al. 2015). Another family of postprocessing 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 signaltonoiseratio (S/N) maps. In contrast with the forward modelbased algorithms, which provide a S/N map as a byproduct of the model, in the case of reference PSF subtraction the S/N map is usually generated by using the medianaveraged residual frames to estimate the annuluswise 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 PSFsubtracted 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 postprocessing. Instead of averaging the set of derotated 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 regimeswitching algorithm to classify the pixel values into two categories, regrouping either the planetary signals or the quasistatic speckles. The probability associated with the planetary regime then allows the creation of a detection map. The algorithm derives fromthe Markov regimeswitching 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 postprocessing methods. The cubes of residuals obtained from the different postprocessing 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 regimeswitching 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 stateoftheart ADIbased postprocessing techniques. Finally, Sect. 5 concludes on this work.
2 Regimeswitching model
The proposed detection algorithm derives from the Markovswitching 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 regimeswitching model (RSM). This approach is one of the most popular nonlinear 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 regimeswitching 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 predefined 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 derotated cube of residuals obtained after the PSF subtraction and derotation steps of the ADI sequence postprocessing. 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 derotated 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 offaxis PSF^{1} or as a forward model of the offaxis PSF after the subtraction. We consider in this paper the measured offaxis PSF for simplicity, although the algorithm may be easily adapted to a forward modelled offaxis PSF.
The RSM we propose here is a modified version of the original Markovswitching model, in which only one parameter is determined via a maximum loglikelihood 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 annuluswise. For each annulus a, a specific residuals time series is built by vectorizing that part of the cube of residuals, indexed by i_{a} the flattened pixel number. The length of the time series depends on the number of pixels in the considered annulus L_{a} but also onthe number of frames in the original derotated cube of residuals T. We indeed take advantage of all the individual frames contained in the derotated 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 i_{a} ∈ {1, …, T × L_{a}}. The first subscript of X indicates the selected frame in the derotated 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 i_{a} 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 regimeswitching model during T steps, allowing theprobability of being in the planetary regime to build up thanks to the shortterm 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 quasistatic 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 (offaxis PSF). The PSF being twodimensional, 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 forwardmodelled offaxis PSF to take into account the signal selfsubtraction, 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 quasistatic speckle residuals, and their time and space varying part characterised by the quasistatic 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 offaxis 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 twostate Markov chain allowing shortterm memory. This implies that we only consider the state in which the system was at index i_{a} − 1 to define the probability of being in agiven state for the current index i_{a}. 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 probabilityweighted sum of the values generated by the equation describing each regime.
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. 
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. 
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 quasistatic 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 i_{a} 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 derotated 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 twostate Markov chain, the computation of necessitates the estimation of (i) the probability of observing the system in the state q at step i_{a} − 1, (ii) the transition probability p_{q,s} from state q to state p and (iii) the likelihood of observing in state s at step i_{a}, which we note . The probability of being in a state s at index i_{a} can be computed as the normalised likelihood of being in state s at index i_{a} multiplied by the probability of having been in either of the two states at index i_{a} −1 and by the transition probability p_{q,s}, which accounts for the shortterm 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 i_{a} 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 i_{a} − 1 via the sum over q. The function , which represents the numerator summed over the two possible states taken at index i_{a}, ensures that the sum of the probability equals one for every index i_{a}.
Description of the mathematical notations for the variables used in the RSM detection map computation.
2.4 Transition probabilities estimation
For our tworegime model, the transition probability p_{q,s} regroups the probabilities of staying in either regime along with the probabilities of switching to the other regime. The estimation of p_{q,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 L_{a} and the number of frames T, the parametrisation of p_{q,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 quasistatic 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 quasistatic speckle residuals. Considering the small transition probabilities p_{0,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 postprocessing algorithms provide different noise distributions for different separations. The Gaussian distribution allows us to construct a likelihood function for state s at index i_{a} 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(S_{t} = q∣P, μ, β, σ) equal to the unconditional probability ξ_{q,0} = P(S_{t} = 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 p_{q,s}, I_{2×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 quasistatic speckles residuals and its first two moments, the planetary signal model P, the intensity parameter β, and the transition probability p_{q,r}. The transition probability p_{q,r} is already defined in Sect. 2.4. We therefore consider the remaining three model parameters.
3.1 Computation of derotated cubes of residuals
The first step to create a RSM detection map is the production of the derotated cubes of residuals for the selected ADIbased postprocessing techniques feeding our regimeswitching 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 postprocessing 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 V_{k}, 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 nonnegativity condition. This method consists in the decomposition of a matrix into two factors of nonnegative 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, M −WH.
Finally, the LLSG estimation is based on the decomposition of the cube intensities in three separate components: L, a lowrank 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 postprocessing 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 postprocessing 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 regimeswitching 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 pixelwise 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 loglikelihood 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.
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. 
3.4 Planetary signal model
Using a forward model for the planetary signal would allow us to take into account the distortions (such as selfsubtraction) created by ADIbased postprocessing treatment when estimating cubes of residuals. Although a forwardmodelled 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 offaxis PSF. We therefore decided to only consider measured offaxis 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 derotate all the resulting frames;
 2.
define the separation to the star for the first and last annuli, respectively a_{ini} = FWHM∕2 + 1 and a_{fin} = (f_{size} − FWHM)∕2 with f_{size} 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 i_{a} for the set of tested δ;
 6.
include the probability of planetary signal providing the maximum likelihood in a threedimensional matrix ;
 7.
repeat steps 2–6 for the next annulus (a + 1) until a_{fin} 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 arcsecond.
The second data set is an ADI sequence on 51 Eridani produced by the SPHEREIRDIS 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 preprocessed 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 postprocessing algorithms. The postprocessing 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 annuluswise, with each annulus being divided into four segments in the case of LLSG. Other parametrisations are possible as the proposed approach works with any derotated cube of residuals. The three cubes of residuals obtained with the selected postprocessing 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 SPHEREIRDIS 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. Fiftyone 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 winddriven 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 ADIbased postprocessing 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.
Fig. 4 Probability map obtained for the SPHEREIRDIS 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. 
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 stateoftheart 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 postprocessing 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 offaxis 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 derotated mediancombined individual frame by estimating the S/N for every pixel. This estimation is done annuluswise 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 postprocessing 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 SPHEREIRDIS data set. As regards the RSM, the mean and variance of the residuals distribution are again estimated annuluswise. 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 loglikelihood 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 postprocessing 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.
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. 
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 paper^{2}. The function allows the user (i) to select one of the two distributions, (ii) to automatically select the best distribution based on a bestfit 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 bestfit approach.
4.3.2 Comparison with S/Nbased detection
We now address the question of the performance of our algorithm compared to the three postprocessing 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 regimeswitching feature. This allows our model to take advantage of the strength of the different postprocessing methods used to produce the cubes of residuals. As speckles are not treated equally by these postprocessing 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 p_{q,s} and on the probabilitiesat step i_{a} − 1 (see Eq. (3)) partly mitigates the effect of speckles on the detection map. Outliers caused by quasistatic 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.
Fig. 6 Detection map obtained after injecting three fake companions in the SPHEREIRDIS 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. 
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 ADIbased postprocessing techniques. The RSM algorithm can be associated with any ADIbased postprocessing 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 quasistatic speckles to be treated more effectively when using several cubes of residuals provided by different postprocessing 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 loglikelihood 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 stateoftheart 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.
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. 
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. 
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
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. 
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. 
References
 Absil, O., Milli, J., Mawet, D., et al. 2013, A&A, 559, L12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Beuzit, J.L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bottom, M., Ruane, G., & Mawet, D. 2017, Res. Notes AAS, 1, 30 [NASA ADS] [CrossRef] [Google Scholar]
 Bowler, B. P. 2016, PASP, 128, 102001 [NASA ADS] [CrossRef] [Google Scholar]
 Cantalloube, F., Mouillet, D., Mugnier, L. M., et al. 2015, A&A, 582, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cosslett, S. R., & Lee, L.F. 1985, J. Econom., 27, 79 [Google Scholar]
 Delorme, P., Meunier, N., Albert, D., et al. 2017, SF2A2017: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, 347 [Google Scholar]
 Goldfeld, S. M., & Quandt, R. E. 1973, J. Econom., 1, 3 [Google Scholar]
 Gomez Gonzalez, C., Absil, O., Absil, P.A., et al. 2016, A&A, 589, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gomez Gonzalez, C. A., Wertz, O., Absil, O., et al. 2017, AJ, 154, 7 [NASA ADS] [CrossRef] [Google Scholar]
 Gomez Gonzalez, C. A., Absil, O., & Van Droogenbroeck M. 2018, A&A, 613, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hamilton, J. D. 1988, J. Econ. Dyn. Control, 12, 385 [Google Scholar]
 Hamilton, J. D. 1994, Time Series Analysis (Princeton: Princeton University Press) [Google Scholar]
 Lafreniere, D., Marois, C., Doyon, R., Nadeau, D., & Artigau, E. 2007, ApJ, 660, 1 [Google Scholar]
 Lagrange, A.M., Bonnefoy, M., Chauvin, G., et al. 2010, Science, 329, 57 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Lozi, J., Guyon, O., Jovanovic, N., et al. 2018, Adaptive Optics Systems VI (Bellingham, USA: SPIE Press) [Google Scholar]
 Macintosh, B. A., Graham, J. R., & Palmer, D. W. 2008, SPIE Conf. Ser., 7015, 701518 [Google Scholar]
 Maire, A.L., Rodet, L., Cantalloube, F., et al. 2019, A&A, 624, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Marois, C., Lafreniere, D., Doyon, R., Macintosh, B., & Nadeau, D. 2006, ApJ, 641, 1 [Google Scholar]
 Mawet, D., Milli, J., Wahhaj, Z., et al. 2014, ApJ, 792, 97 [NASA ADS] [CrossRef] [Google Scholar]
 Pairet, B., Cantalloube, F., Gomez Gonzalez, C., Absil, O., & Jacques, L. 2019, MNRAS, 487, 2262 [NASA ADS] [CrossRef] [Google Scholar]
 Pueyo, L. 2016, ApJ, 824, 117 [NASA ADS] [CrossRef] [Google Scholar]
 Ren, B., Pueyo, L., Zhu, G. B., Debes, J., & Duchêne, G. 2018, ApJ, 852 [Google Scholar]
 Ruffio, J.B., Macintosh, B., Wang, J. J., & Pueyo, L. 2017, ApJ, 842, 1 [Google Scholar]
 Samland, M., Mollière, P., Bonnefoy, M., et al. 2017, A&A, 603, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Soummer, R., Pueyo, L., & Larkin, J. 2012, ApJ, 755, L2 [Google Scholar]
 Vigan, A., Gry, C., Salter, G., et al. 2015, MNRAS, 454, 129 [NASA ADS] [CrossRef] [Google Scholar]
For coronagraphic imaging, an offaxis noncoronagraphic 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 modelbased algorithm. For noncoronagraphic imaging, this reference PSF is the unsaturated exposure.
The RSM detection map python package is available on GitHub: https://github.com/chdahlqvist/RSMmap
All Tables
Description of the mathematical notations for the variables used in the RSM detection map computation.
All Figures
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. 

In the text 
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. 

In the text 
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. 

In the text 
Fig. 4 Probability map obtained for the SPHEREIRDIS 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. 

In the text 
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. 

In the text 
Fig. 6 Detection map obtained after injecting three fake companions in the SPHEREIRDIS 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. 

In the text 
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. 

In the text 
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. 

In the text 
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. 

In the text 
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. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.