# deep PACO: Combining statistical models with deep learning for exoplanet
detection and characterization in direct imaging at high contrast

Olivier Flasseur1,2,5, Théo Bodrito2, Julien Mairal3, Jean Ponce2,4, Maud
Langlois5, and Anne-Marie Lagrange1,6
1Laboratoire d’Études Spatiales et d’Instrumentation en Astrophysique,
Observatoire de Paris, Univ. PSL, Sorbonne Univ., Univ. Paris Diderot, France
2Département d’Informatique de l’École Normale Supérieure (ENS-PSL, CNRS,
Inria), France
3Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK, 38000 Grenoble, France
4Courant Institute of Mathematical Sciences, Center for Data Science, New York
Univ., USA
5Centre de Recherche Astrophysique de Lyon, CNRS, Univ. de Lyon, Univ. Claude
Bernard Lyon 1, ENS de Lyon, France
6Univ. Grenoble Alpes, Institut de Planétologie et d’Astrophysique de
Grenoble, France O. Flasseur is currently with CRAL (5), and was previously
with LESIA (1) and Inria (2) where most of this work was performed. E-mail:
<EMAIL_ADDRESS>
###### Abstract

Direct imaging is an active research topic in astronomy for the detection and
the characterization of young sub-stellar objects. The very high contrast
between the host star and its companions makes the observations particularly
challenging. In this context, post-processing methods combining several images
recorded with the pupil tracking mode of telescope are needed. In previous
works, we have presented a data-driven algorithm, PACO, capturing locally the
spatial correlations of the data with a multi-variate Gaussian model. PACO
delivers better detection sensitivity and confidence than the standard post-
processing methods of the field. However, there is room for improvement due to
the approximate fidelity of the PACO statistical model to the time evolving
observations. In this paper, we propose to combine the statistical model of
PACO with supervised deep learning. The data are first pre-processed with the
PACO framework to improve the stationarity and the contrast. A convolutional
neural network (CNN) is then trained in a supervised fashion to detect the
residual signature of synthetic sources. Finally, the trained network delivers
a detection map. The photometry of detected sources is estimated by a second
CNN. We apply the proposed approach to several datasets from the VLT/SPHERE
instrument. Our results show that its detection stage performs significantly
better than baseline methods (cADI, PCA), and leads to a contrast improvement
up to half a magnitude compared to PACO. The characterization stage of the
proposed method performs on average on par with or better than the comparative
algorithms (PCA, PACO) for angular separation above 0.5”.

###### keywords:

techniques: high angular resolution – techniques: image processing – methods:
numerical – methods: statistical – methods: data analysis

††pubyear: 2023††pagerange: deep PACO: Combining statistical models with deep
learning for exoplanet detection and characterization in direct imaging at
high contrast–B

##  1 Introduction

High-contrast imaging is an observational method used to study the close
environment of stars (Traub & Oppenheimer, 2010; Bowler, 2016; Pueyo, 2018).
It is particularly adapted to detect young, massive and hot exoplanets (see
e.g. Chauvin et al. (2004); Chauvin et al. (2005); Schneider et al. (2011);
Nielsen et al. (2019)), thus complementing well indirect exoplanet detection
methods such as transit photometry or Doppler spectroscopy (Santos, 2008).
Direct imaging offers other appealing characteristics like the detection of
candidate companions from a few hours of observations and the ability to
characterize them in terms of age, surface gravity, effective temperature and
composition (Allard et al., 2003, 2007), or to predict their evolution
(Burrows et al., 1997; Chabrier et al., 2000). Despite these promises, only a
few dozens exoplanets have been unveiled and characterized since the emergence
of direct imaging in the early 2000s (Marois et al., 2008; Lagrange et al.,
2009; Nielsen et al., 2012; Macintosh et al., 2015; Chauvin et al., 2017;
Keppler et al., 2018). This is mainly due to (i) the relatively low occurrence
of giant exoplanets, (ii) the very high contrast between the host star and the
exoplanets (typically, higher than $10^{5}$ in the infrared), and (iii) the
required high angular resolution (typically, better than a tenth of arcsec).

In this context, cutting-edge ground-based facilities like VLT/SPHERE (Beuzit
et al., 2019), Gemini/GPI (Macintosh et al., 2014), Keck/NIRC2 (Castellá et
al., 2016), Magellan/MagAO (Morzinski et al., 2014) or SUBARU/SCExAO
(Jovanovic et al., 2015) are equipped with an (extreme) adaptive optics system
and a coronagraph to attenuate as much as possible the glare of the star.
Currently, the non-blocked residual starlight contamination and its temporal
fluctuations remain the main limitation. It takes the form of spatially-
correlated speckles that resemble the expected signature of a point-like
source (e.g., an exoplanet, a brown dwarf, a background star). The
observations are also impacted by additional sources of noise (i.e., thermal
background flux, detector readout, photon noise). Together, speckles and noise
form a spatially and non-stationary nuisance component that corrupts the
signals of the sought objects. Off-axis objects can either take the form of
point-like sources or that of spatially extended features like circumstellar
disks. In this paper, we focus on the detection of point-like sources, leaving
the problem of disk reconstruction for future work.

In order to unmix the sought objects from the nuisance component, high-
contrast observations are performed with dedicated strategies like angular
differential imaging (ADI; Marois et al. (2006)), that we consider in this
paper. ADI consists in tracking the observed target over time, with the
telescope derotator tuned to keep the telescope pupil stable while the field
of view rotates around the host star. Consequently, in the resulting 3-D
datasets (2-D + time), the objects of interest follow an apparent motion along
a deterministic circular trajectory centered on the star while the telescope
pupil remains static. With ADI, speckles due to residual starlight aberrations
are quasi-static, i.e., they are strongly correlated across exposures. ADI
also allows to extract the astrometry and the photometry of detected sources.
These estimates can be used to characterize the physical properties of the
detected sources by fitting orbital and atmospheric models (Vigan et al.,
2010; Cheetham et al., 2019; Mesa et al., 2019). In this paper, we address two
tasks: (i) the detection of point-like sources, and (ii) the estimation of
their photometry.

The last cornerstone of high-contrast imaging is data reduction, i.e., the
combination of the recorded images by dedicated post-processing algorithms.
This critical step brings the additional gain in contrast (between one and
three orders of magnitudes for existing methods) needed to detect the faint
signals coming from thermal self-emission of giant exoplanets. The classical
principle is to estimate a reference image (so-called on-axis PSF) of the
nuisance component, that can be subtracted from the data in order to unveil
the sought objects. A simple practical implementation of this general strategy
consists in subtracting the temporal mean or median of the dataset from each
frame of the ADI stack. The residual images are then co-aligned to the true-
North so that the signals of the sought objects are superimposed and can be
combined by temporal stacking. This is the principle of the cADI method
designed to process the first direct imaging observations (Marois et al.,
2006; Lagrange et al., 2009). In the last decade, several more advanced
methods have been developed, see, e.g., Pueyo (2018) for a review. In
particular, TLOCI (Marois et al., 2013; Marois et al., 2014) (or its variants
such as LOCI (Lafrenière et al., 2007), ALOCI (Currie et al., 2012a; Currie et
al., 2012b), MLOCI (Wahhaj et al., 2015)), and KLIP/PCA (Soummer et al., 2012;
Amara & Quanz, 2012) are currently implemented in most reduction pipelines
(Amara & Quanz, 2012; Gonzalez et al., 2017; Galicher et al., 2018). LOCI-
based and PCA-based algorithms are considered as standards to process high-
contrast observations. With the (A, M, T)-LOCI algorithm, the on-axis PSF is
estimated by combining images selected in a library. The combined images are
selected and weighted to minimize the residual noise and maximize the
throughput of point-like sources simultaneously. The PCA algorithm performs a
principal component analysis of the data, and a low-rank estimate of the on-
axis PSF is formed by keeping the first principal components of the
decomposition. In the same vein, the LLSG algorithm (Gonzalez et al., 2016)
decomposes the dataset into low-rank, sparse and Gaussian components. Because
few sources are expected in the field of view, their signatures are mainly
recovered in the sparse component. The RSM algorithm (Dahlqvist et al., 2020,
2021a, 2021b) combines residual images obtained with different post-processing
algorithms (e.g., cADI and PCA) to leverage their specific benefits and
mitigate their respective drawbacks. Since all of these algorithms are based
on the estimation and subtraction of an off-axis PSF, they are facing a common
pitfall: they fail in deriving a statistically grounded detection map,
especially at short angular separations. As a consequence, the identification
of candidate sources partly relies on visual inspection of the detection map.

To circumvent this issue, several works have considered alternative ways to
reduce the data in order to produce more quantitative outputs. The derivation
of a custom signal-to-noise ratio (S/N) through a $t$-test empirically
corrected for the varying number of samples as a function of the angular
separation (Mawet et al., 2014) is a pioneering work in this direction.
Jensen-Clem et al. (2017) recommend the adoption of metrics combining the
achievable contrast with the fraction of detected sources. Instead of
normalizing the combined residual images by the empirical standard-deviation
of the noise (roughly) approximated on annuli, Pairet et al. (2019) propose to
build a detection map directly from the set of residual images by comparing
the variance of samples located on the expected trajectory of putative
sources. Other methods reformulate the detection task as an inverse problem.
Among them, ANDROMEDA (Mugnier et al., 2009; Cantalloube et al., 2015) and
FMMF (Ruffio et al., 2017) build a model of the residual off-axis signal after
subtraction of the estimated on-axis PSF. The PACO algorithm (Flasseur et al.,
2018a, b, c, 2020a, 2020b) builds a more consistent statistical model, self-
calibrated on the data, that accounts for the spatial correlations of the
nuisance component at the scale of small image patches of a few tens of
pixels.

Given the success of data-driven approaches in solving various high-level
imaging tasks (e.g., detection, segmentation, classification) in very diverse
fields (e.g., photography, microscopy, bio-medical imaging, remote sensing),
machine learning and deep learning approaches have also been investigated by
the direct imaging community. Fergus et al. (2014) report one of the first
steps in this direction with a discriminative approach exploiting the specific
structure of high-contrast data. Based on support vector machines, the
underlying model is trained from two-classes samples generated by resorting to
massive injections of fake companions. Gonzalez et al. (2018) formalize the
detection problem as a binary classification task and propose a fully
supervised deep learning approach also trained in a supervised fashion. They
use collections of patches pre-processed by PCA for different numbers of
principal components as input of a random forest or of a convolutional neural
network that decides in favor of the presence or on the absence of a point-
like source in each patch. While demonstrating powerful detection
capabilities, this algorithm showed to be prone to a high level of false
alarms in some cases (Cantalloube et al., 2020). Besides, the tuning of hyper-
parameters remains a critical point making the operating point difficult to
reach. Generative adversarial networks (GANs, Goodfellow et al. (2014)) have
been used to produce multiple realizations of pure nuisance component as an
alternative way to generate a large basis of labeled samples used to train a
deep discriminative model (Yip et al., 2019). Recent works (Samland et al.,
2021; Gebhard et al., 2022) recast the unmixing problem as a regularized
regression task. The underlying linear model is in charge of explaining the
evolution of the nuisance component (and possibly, the evolution of the source
signals) in a time series extracted at a given pixel location from temporal
series of reference selected to be signal free and causally independent from
the putative source signals.

Intensive testing of PACO, both on public (Cantalloube et al., 2020) as well
as on consortia data challenges, and more recently on a sub-sample of about 75
observations from the SPHERE’s SHINE survey (Chomez et al., 2022) shows that
PACO is one of the algorithms of choice to process high-contrast observations.
In particular, the latter work (Chomez et al., 2022) demonstrates the expected
gain in terms of achievable contrast (up to $10^{-7}$ at a few arcsec), and in
terms of the underlying exoplanet population (with a mass up to
$5\,\text{M}_{\text{Jup}}$ at 5 AU for stars at about 60 parsec away). Thanks
to its unique data-driven modeling of the nuisance component, accounting for
its non-stationary spatial correlations, PACO is especially well suited to
process observations in which the typical spatial extent of speckles lies in a
patch of a few tens of pixels. However, the statistical model of PACO is
approximate in case of spatial correlations spread over a patch of a few tens
of pixels (e.g., for background-limited observations and/or in case of
unstable observing conditions). This motivates the path we follow in this
work: we propose a new detection algorithm, named deep PACO, that combines the
statistical model of PACO with a supervised deep learning framework. The
statistical model of PACO is used to improve the stationarity and the contrast
of the data in a pre-processing step, while deep learning is in charge of
correcting for the (putative) approximate fidelity of the statistical model of
PACO to the reality of the observations. To do so, the data are centered and
whitened locally using the PACO framework, and a CNN is trained in a
supervised fashion to detect the residual signature of synthetic sources from
pre-processed science data. The network is trained from scratch using full
frame samples generated with a custom data augmentation strategy allowing to
build a large training set from a single ADI dataset. Finally, the underlying
discriminative model is applied to the pre-processed observations and delivers
a detection map. Detected sources are then photometrically characterized by a
second deep neural network, also trained from scratch using patch samples
generated with a dedicated data augmentation step. On this latter part, while
the recent astronomy literature reports several works for photometry
estimation through deep learning models (see e.g., Boucaud et al. (2020);
Cabayol et al. (2021); Huertas-Company & Lanusse (2022) for galaxie’s
photometry or red-shift estimation), to the best of your knowledge this is the
first time that deep learning is employed to estimate the flux of detected
sources in direct imaging at high-contrast.

This paper111A preliminary version of this work was presented in the form of a
conference contribution in Flasseur et al. (2022). The present manuscript
contains a significant amount of additional methodological developments,
technical improvements, and experiments. is organized as follow. Section 2
presents the main ingredients of the detection part of the proposed algorithm.
Section 3 focuses on the characterization stage of the proposed method.
Section 4 evaluates the detection and characterization performance on several
high-contrast observations from the IRDIS imager (Dohlen et al., 2008) of the
VLT/SPHERE instrument (Beuzit et al., 2019). Finally, Sect. 5 draws the paper
conclusions and gives future research prospects.

Throughout the text, the reader can refer to Table 1, Fig. 1, and Table 11,
summarizing respectively the main notations, the processing pipeline of the
proposed approach, and the main setting of the different parameters.

Table 1: Summary of the main notations. Not. | Range | Definition
---|---|---
$\blacktriangleright$ constants
$N$ | $\mathbb{N}$ | number of pixels in a frame
$M$ | $\mathbb{N}$ | number of pixels in a detection map
$T$ | $\mathbb{N}$ | number of temporal frames
$K$ | $\mathbb{N}$ | number of pixels in a patch (pre-processing)
$J$ | $\mathbb{N}$ | number of pixels in a patch (characterization)
$Q$ | $\mathbb{N}$ | number of samples involved in shrinkage
$P$ | $\mathbb{N}$ | total number of training sources
$S$ | $\mathbb{N}$ | total number of training sets
$\blacktriangleright$ indexes
$n$ | $\llbracket 1;N\rrbracket$ | pixel index
$t$ | $\llbracket 1;T\rrbracket$ | temporal index
$p$ | $\llbracket 1;P\rrbracket$ | source index
$s$ | $\llbracket 1;S\rrbracket$ | training set index
$\phi$ | $\mathbb{R}_{+}^{2}$ | 2-D (sub-pixel) angular location of a source
$\blacktriangleright$ data quantities
$\bm{r}$ | $\mathbb{R}^{N\times T}$ | observations
$\bm{f}$ | $\mathbb{R}^{N\times T}$ | nuisance component
$\bm{h}$ | $\mathbb{R}^{N}$ | off-axis PSF
$\widetilde{\bm{r}}$ | $\mathbb{R}^{N\times T}$ | pre-processed observations without injections
$\widebar{\bm{r}}$ | $\mathbb{R}^{N\times T}$ | observations with injections
$\widecheck{\bm{r}}$ | $\mathbb{R}^{N\times T}$ | pre-processed observations with injections
$\mathchoice{\mathop{\vbox{\halign{#\cr\kern 0.80002pt$\hss\leavevmode\resizebox{3.8315pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\displaystyle\bm{r}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern 0.80002pt$\hss\leavevmode\resizebox{3.8315pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\textstyle\bm{r}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern 0.80002pt$\hss\leavevmode\resizebox{2.68205pt}{2.72223pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptstyle\bm{r}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern 0.80002pt$\hss\leavevmode\resizebox{1.91574pt}{1.94444pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptscriptstyle\bm{r}\hss$\crcr}}}\limits}$ | $\mathbb{R}^{N\times T}$ | input (images) of the CNN (detection)
$\mathchoice{\mathop{\vbox{\halign{#\cr\kern 0.80002pt$\hss\leavevmode\resizebox{4.02501pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\displaystyle\bm{p}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern 0.80002pt$\hss\leavevmode\resizebox{4.02501pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\textstyle\bm{p}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern 0.80002pt$\hss\leavevmode\resizebox{2.8175pt}{2.72223pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptstyle\bm{p}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern 0.80002pt$\hss\leavevmode\resizebox{2.0125pt}{1.94444pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptscriptstyle\bm{p}\hss$\crcr}}}\limits}$ | $\mathbb{R}^{J}$ | input (patch) of the CNN (characterization)
$\blacktriangleright$ operators
$\mathbf{E}$ | $\mathbb{R}^{N\times\,.}\mapsto\mathbb{R}^{K\times\,.}$ | patch extractor (pre-processing)
$\mathbf{W}$ | $\mathbb{R}^{K\times\,.}\mapsto\mathbb{R}^{K\times\,.}$ | centering and whitening (pre-processing)
$\mathbf{D}$ | $\mathbb{R}^{\,.\times T}\mapsto\mathbb{R}^{\,.\times T}$ | frame derotator by parallactic angles
$\blacktriangleright$ losses and metrics
$\mathcal{L}_{\text{detect.}}$ | $\left[0;1\right]^{2\times M}\mapsto\mathbb{R}_{+}$ | Dice2 score (detection loss)
$\mathcal{L}_{\text{carac.}}$ | $\mathbb{R}_{+}^{2}\mapsto\mathbb{R}_{+}$ | absolute relative error (characterization loss)
TPR | $\left[0;1\right]$ | true positive rate (detection metric)
FDR | $\left[0;1\right]$ | false discovery rate (detection metric)
F1R | $\left[0;1\right]$ | harmonic mean of TPR and FDR (detection metric)
$\blacktriangleright$ estimated quantities
$\widehat{\bm{y}}$ | $\left[0;1\right]^{M}$ | detection map
$\widehat{\alpha}$ | $\mathbb{R}_{+}$ | photometry (in source to star contrast)
$\widehat{\bm{m}}$ | $\mathbb{R}^{N}$ | temporal mean
$\widehat{\rho}$ | $\left[0;1\right]$ | shrinkage factor
$\widehat{\mathbf{S}}$ | $\mathbb{R}^{K\times K}$ | empirical spatial covariance
$\widehat{\mathbf{C}}$ | $\mathbb{R}^{K\times K}$ | shrunk spatial covariance
$\widehat{\mathbf{L}}$ | $\mathbb{R}^{K\times K}$ | Cholesky’s factorization of $\widehat{\mathbf{C}}$
Figure 1: (a) Typical observations $\bm{r}$ from the VLT/SPHERE-IRDIS
instrument conducted in pupil tracking mode (i.e., with the ADI technique).
Zooms around two known background faint sources are displayed in the red
circles, a zoom on the nuisance component $\bm{f}$ near the star is
highlighted in the green circle, and a view of the sought exoplanetary signal,
taking the form of the off-axis PSF $\bm{h}$, is shown in the blue circle. Two
spatio-temporal slice cuts along the black solid and dashed lines are shown on
the right, as an illustration of the spatial non-stationarity of the
correlations of the nuisance. (b) Illustration of the main operations
performed during the pre-processing of the proposed algorithm by statistical
learning. Examples of estimated spatial covariance matrices are displayed in
squares for four regions of interest pickled at about $0.5$ (pink), 1.0
(green), 1.5 (purple), and 2.0 (orange) arcsec. (c) Schematic representation
of the main operations performed during the detection and characterization
steps of the proposed algorithm by supervised deep learning. The first line
displays common operations performed for both the detection and the
characterization. The left (respectively, the right) part of the second line
is for operations applied solely during the detection (respectively, the
characterization) step. Dataset: HIP 72192 (2015-05-10), see Sect. 4 for the
recording logs.

##  2 Source detection algorithm

###  2.1 Stastical model of the non-stationary patch covariances

Section 2.1.1 recalls for completeness the main ingredients of the statistical
model embedded in the PACO algorithm for ADI observations (Flasseur et al.,
2018b, 2020a). Section 2.1.2 describes how to use this model to attenuate the
strong and spatially non-stationary correlations of the data as well as
improve the contrast. Ablation studies have shown that this pre-processing
step is of primary importance: the training step fails to converge in its
absence due to the high non-stationarity of the nuisance component and the
large fluctuations it displays near the star, see Fig. 1 (a)-(b). This effect
is due to the fact that standard deep learning architectures assume some
degree of invariance, in particular that the data are normalized in a specific
range and are corrupted by a stationary noise.

####  2.1.1 Statistical model of the nuisance component

A dataset $\bm{r}$ in $\mathbb{R}^{N\times T}$ recorded with the ADI technique
is formed by $N$-pixel images captured at different times $t$ in $\llbracket
1;T\rrbracket$. The direct model for the observed intensity is:

$\bm{r}=\bm{f}+\sum\limits_{p=1}^{P}\alpha_{p}\,\bm{h}(\phi_{p})\,,$ (1)

where $\bm{f}$ in $\mathbb{R}^{N\times T}$ is the nuisance component, and
$\bm{h}\left(\phi_{p}\right)$ in $\mathbb{R}^{N\times T}$ stands for the
contribution of a point source $p\in\llbracket 1;P\rrbracket$ with a contrast
$\alpha_{p}$ in $\mathbb{R}_{+}$ that is assumed constant during the few hours
of the total observations. The contribution of a source $p$ takes the form of
the off-axis PSF centered at location $\mathcal{F}_{t}(\phi_{p})$ in the
$t$-th image, where $\phi_{p}$ is its initial location on an image at a
reference time $t_{\text{ref}}$ (e.g., $t_{\text{ref}}=t_{1}$). The function
$\mathcal{F}_{t}$ is a geometric transform (typically in ADI, a circular
translation with respect to the star located at the center of the images)
modeling the apparent motion of the field of view between the observation
configurations at time $t_{\text{ref}}$ and time $t$. The mapping
$\mathcal{F}_{t}$ is deterministic since it depends solely on the measured
parallactic angles. Given that very few sources are expected in the field of
view, we assume that the measured intensity is the superimposition of the
nuisance component and at most one unresolved point-like source $p$ at each
pixel location $n$, i.e., multiple sources do not overlap.

In previous works on the PACO algorithm (Flasseur et al., 2018a, b, c), we
have proposed to describe the random fluctuations of the nuisance component
$\bm{f}$ by a statistical model whose parameters are estimated in a data-
driven fashion. We recall hereafter the main ingredients of this statistical
model.

Given the spatial non-stationarity of the nuisance component, the model is
built locally at the scale of a patch with an area of a few tens of pixels. It
models the distribution of $T$ patches222For the convenience, we use in the
equations a vectorized version of 2-D patches.
$\bm{f}_{n}=\\{\mathbf{E}_{n,t}\,\bm{f}\\}_{t=1:T}$ in $\mathbb{R}^{K\times
T}$ extracted around pixel $n$ ($\mathbf{E}_{n,t}$ denotes the $K$-pixel patch
extraction operator at location $n$ and time $t$) with a multi-variate
Gaussian $\mathcal{N}(\bm{m}_{n},\mathbf{C}_{n})$. The covariance matrix
$\mathbf{C}_{n}$ is non-diagonal, i.e., it accounts for the local correlations
of $\bm{f}$. The sample estimators
$\\{\widehat{\bm{m}}_{n};\widehat{\mathbf{S}}_{n}\\}$ of the local mean and
covariances coming from the maximum likelihood are the following:

$\begin{cases}\widehat{\bm{m}}_{n}=\frac{1}{T}\sum\limits_{t=1}^{T}\mathbf{E}_{n,t}\,\bm{r}\,\in\mathbb{R}^{K}\,,\vspace{0.5mm}\\\
\widehat{\mathbf{S}}_{n}=\frac{1}{T}\sum\limits_{t=1}^{T}(\mathbf{E}_{n,t}\,\bm{r}-\widehat{\bm{m}}_{n})(\mathbf{E}_{n,t}\,\bm{r}-\widehat{\bm{m}}_{n})^{\top}\,\in\mathbb{R}^{K\times
K}\,.\end{cases}$ (2)

Since the number $T$ of temporal frames is typically equivalent or smaller
than the number $K$ of pixels in a patch333The number of pixels in a patch is
determined in a data-driven fashion, as described in Flasseur et al. (2018b)
for the PACO algorithm. It corresponds to twice the full width at half maximum
(FWHM) of the measured off-axis PSF at the given wavelength. In practice, this
empirical rule typically leads to $K$ in $\llbracket 7^{2};12^{2}\rrbracket$
pixels for the VLT/SPHERE instrument operating at a wavelength
$\lambda\in\left[0.9;2.2\right]\,\mathrm{\SIUnitSymbolMicro}\mathrm{m}$., the
sample covariance $\widehat{\mathbf{S}}_{n}$ is very noisy and may even be
rank deficient. A form of regularization must be enforced to stabilize the
estimate and allow the inversion of the covariance matrix involved in the data
whitening step (see Sect. 2.1.2). We use a shrinkage estimator (Ledoit & Wolf,
2004; Chen et al., 2010) formed by a convex combination between the low
bias/high variance estimator $\widehat{\mathbf{S}}_{n}$ and a high bias/low
variance estimator $\widehat{\mathbf{F}}_{n}$:

$\widehat{\mathbf{C}}_{n}=(1-\widehat{\rho}_{n})\,\widehat{\mathbf{S}}_{n}+\widehat{\rho}_{n}\,\widehat{\mathbf{F}}_{n}\,,$
(3)

where $\widehat{\mathbf{F}}_{n}$ is a diagonal matrix encoding the sample
variances:

$\left[\widehat{\mathbf{F}}_{n}\right]_{kk^{\prime}}=\begin{cases}\left[\widehat{\mathbf{S}}_{n}\right]_{kk^{\prime}}~{}~{}~{}~{}~{}&\text{if}~{}~{}~{}~{}~{}k=k^{\prime}\\\
0~{}~{}~{}~{}~{}&\text{if}~{}~{}~{}~{}~{}k\neq k^{\prime}\,.\end{cases}$ (4)

The hyper-parameter $\widehat{\rho}_{n}$ plays a central role since it governs
a bias-variance trade-off. In our previous works (Flasseur et al., 2018b,
2021), we have derived its closed-form expression, which is an extension of
the results of Chen et al. (2010) in the case of a non-constant shrinkage
matrix $\widehat{\mathbf{F}}_{n}$:

$\widehat{\rho}_{n}=\frac{\operatorname{tr}\big{(}\widehat{\mathbf{S}}_{n}^{2}\big{)}+\operatorname{tr}^{2}\big{(}\widehat{\mathbf{S}}_{n}\big{)}-2\sum_{k=1}^{K}\big{[}\widehat{\mathbf{S}}_{n}\big{]}_{kk}^{2}}{(Q+1)\,\Big{(}\operatorname{tr}\big{(}\widehat{\mathbf{S}}_{n}^{2}\big{)}-\sum_{k=1}^{K}\big{[}\widehat{\mathbf{S}}_{n}\big{]}_{kk}^{2}\big{)}}\,,$
(5)

where $Q$ is the number of non-null patches involved in the computation of
$\widehat{\mathbf{S}}_{n}$. Here, $Q$ is equal to $T$ everywhere.

In Appendix A.1, we discuss a refinement of this statistical model to account
for the temporal fluctuations of the observations. It leads to a slight
improvement in terms of detection performance at the cost of an increase of
the computational burden by a factor 10 to 30. For these reasons, it is not
applied by default in the following. We recommend to use it, in a second step,
to refine the analysis of ambiguous candidate detections found by the proposed
method embedding a multi-variate Gaussian model.

####  2.1.2 Centering and local whitening of the observations

We consider a set of locations $\mathbb{P}$ where the statistics of the
nuisance component should be computed. The cardinal of $\mathbb{P}$ depends
solely on the patch shape and on the patch stride used to cover the whole
field of view. By default, we consider square patches of $K$ pixels area. For
a given patch stride, we first define the number of (centered and whitened)
patches averaged at each location $n^{\prime}$ of the field of view:

$q_{n^{\prime}}=\sum\limits_{n\in\mathbb{P}}\delta_{\mathbf{1}_{n}\,\mathbf{1}_{n^{\prime}}^{\top}\,\neq\,0}\,,\forall
n^{\prime}\in\llbracket 1;N\rrbracket\,,$ (6)

where $\delta$ is the indicator function (i.e., $\delta_{x=y}$ is either equal
to 1 if the condition $x=y$ is met, and equal to 0 otherwise), and
$\mathbf{1}_{n}\in\mathbb{R}^{K}$ (respectively,
$\mathbf{1}_{n^{\prime}}\in\mathbb{R}^{K}$) is the vectorization of a 1-valued
patch centered at location $n$ (respectively, $n^{\prime}$). The pre-processed
images $\widetilde{\bm{r}}$ in $\mathbb{R}^{N\times T}$, after centering and
whitening, are obtained by:

$\displaystyle\widetilde{r}_{n^{\prime}}$
$\displaystyle=\frac{\left[\sum\limits_{n\in\mathbb{P}}\mathbf{E}_{n}^{\top}\,\mathbf{W}_{n}\,\mathbf{E}_{n}\,\bm{r}\right]_{n^{\prime}}}{q_{n^{\prime}}}\,,$
$\displaystyle=\frac{\left[\sum\limits_{n\in\mathbb{P}}\mathbf{E}_{n}^{\top}\,\widehat{\mathbf{L}}_{n}^{\top}\,(\bm{r}_{n}-\widehat{\bm{m}}_{n})\right]_{n^{\prime}}}{q_{n^{\prime}}}\,,\forall
n^{\prime}\in\llbracket 1;N\rrbracket\,,$ (7)

where $\mathbf{W}_{n}$ is an operator performing centering and whitening of
the collection of patches $\bm{r}_{n}$ $\mathbb{R}^{K\times T}$ at location
$n$, such that $\widehat{\mathbf{L}}_{n}$ is the Cholesky’s factorization of
$\widehat{\mathbf{C}}_{n}^{-1}$ (i.e.,
$\widehat{\mathbf{L}}_{n}\,\widehat{\mathbf{L}}_{n}^{\top}=\widehat{\mathbf{C}}_{n}^{-1}$).
In the specific case of non-overlapping square patches of $K$ pixels,
$\text{card}(\mathbb{P})=\lfloor N/K\rceil$, and Eqs. (6)-(7) simplify as:

$\begin{cases}q_{n}^{\prime}=1\,,\forall n^{\prime}\in\llbracket
1;N\rrbracket\,,\\\
\widetilde{\bm{r}}_{n}=\mathbf{W}_{n}\bm{r}_{n}=\widehat{\mathbf{L}}_{n}^{\top}\,\left(\bm{r}_{n}-\widehat{\bm{m}}_{n}\right)\,,\forall
n\in\mathbb{P}\,.\end{cases}$ (8)

In Appendix A.2, we discuss a refinement of the intermediate quantity
$\widetilde{\bm{r}}$ produced by the pre-processing step under the same
statistical model as described in Sect. 2.1.1 or in Appendix A.1. It leads to
a better stability and robustness of the method for datasets recorded under
bad observing conditions444Since the notion of observing conditions is
relative and can be characterized by several metrics (e.g., Strehl ratio,
coherence time, air mass, etc.), we did not find quantitative values for these
measures indicating the strict use of the variant version described in
Appendix A.2. This large effort would require the processing of hundreds of
datasets, which is left for future work. Qualitatively, we observed that
datasets impacted by a bright wind-driven halo, displaying the same apparent
motion than the objects of interest, are more subject to lead to an increased
false alarm rate without the pre-processing described in Appendix A.2..
However, the computational burden of this variant is increased by a factor
$2\times K$, typically lying in $\llbracket 100;250\rrbracket$ for the
VLT/SPHERE instrument. For these reasons, this variant is not applied by
default in the following. We recommend to apply it, in a second step, when the
training, validation, or inference results are clearly impacted by a
significant number of false alarms, much higher than expected, which is an
unambiguous sign of the limited fidelity to the observations of the pre-
processing procedure used by default.

###  2.2 Exoplanet detection by supervised deep learning

We formalize the detection problem as a supervised pixel-wise classification
task: starting from a temporal series of pre-processed images including
synthetic sources, the goal is to infer a detection map $\widehat{\bm{y}}$ in
$\left[0;1\right]^{M}$, where each pixel value represents a score between 0
and 1 such that a high (respectively, a low) score values the presence
(respectively, the absence) of a source centered at that location.
Interpreting this score as a true probability of presence of a source requires
a control of the uncertainties with dedicated methods (see e.g., Gawlikowski
et al. (2021); Hüllermeier & Waegeman (2021) for recent review papers) that is
left for future work. For this reason, in the following, we refer to this
score as a pseudo-probability. Besides, the $M$-pixel detection map is larger
than a $N$-pixel single temporal frame of the pre-processed data. Due to the
apparent rotation induced by ADI, it is theoretically possible to detect a
source lying in the sensor field of view at a single date. By combining all
measurements, an area larger than $N$ pixels can be scrutinized (and its
spatial extent depends solely on the maximum parallactic rotation).

Section 2.2.1 details the construction process of the training set, Sect.
2.2.2 describes the selected model architecture, and Sect. 2.2.3 discusses the
metrics we consider to evaluate the performance of the proposed method.

####  2.2.1 Construction of the training basis

In high-contrast imaging, obtaining real ground-truth data is a twofold
challenge.

First, the overall number of positive samples is limited as relatively few
point-like sources have been confirmed to date. Second, negative samples are
hard to define since some undiscovered sources might be present in the
observed data. To overcome these difficulties, we adopt the following training
strategy: the training set consists of $S$ pairs
$\\{{\mathchoice{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{3.8315pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\displaystyle\bm{r}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{3.8315pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\textstyle\bm{r}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{2.68205pt}{2.72223pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptstyle\bm{r}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{1.91574pt}{1.94444pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptscriptstyle\bm{r}\hss$\crcr}}}\limits}^{[s]};\bm{y}^{[s]}}\\}_{s=1:S}$
of samples resulting from the massive injection of synthetic point-like
sources. In this framework, $\mathchoice{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{3.8315pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\displaystyle\bm{r}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{3.8315pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\textstyle\bm{r}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{2.68205pt}{2.72223pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptstyle\bm{r}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{1.91574pt}{1.94444pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptscriptstyle\bm{r}\hss$\crcr}}}\limits}^{[s]}\in\mathbb{R}^{N\times
T}$ represents observations, with injected synthetic sources, that have been
pre-processed. The quantity $\bm{y}^{[s]}\in\llbracket 0,1\rrbracket^{M}$ is
the ground-truth map pointing the injection locations. The implemented
simulation process is quite realistic since the injected source signature
corresponds to the off-axis PSF of the target star usually measured just
before or just after the main sequence of observations by decentering the
coronagraph.

Second, the nuisance component varies drastically from one observation to the
other, as it is highly dependent on the observing conditions, the magnitude of
the star, and the instrument settings. As a consequence, we follow an
observation-dependent approach, and train a different model on each
observation. It means that the model parameters (except algorithmic and
optimization hyper-parameters, see Sect. 2.2.4) are optimized from scratch for
each dataset.

This setup implies the design of a custom data-augmentation strategy (i) to
prevent over-fitting of the model that is trained from a unique temporal
series of images, and (ii) to account for our lack of knowledge about real
sources –unknown at training time but that we aim to detect at test time–. To
circumvent these issues, we apply a random permutation of the $T$ images
forming the observations $\bm{r}$ for each new training sample $s\in\llbracket
1;S\rrbracket$. This operation allows us (i) to create artificially different
training sets, and (ii) to break the temporal consistency of (known and
unknown) real sources. Synthetic sources are then injected in the temporally
permuted data using the parallactic angles and the best fit of the off-axis
PSF by a Gaussian and an Airy pattern. Besides, the off-axis PSF is assumed to
be time-invariant. We have checked numerically that this assumption is
reasonable for our classification task. Similarly, assuming a slightly
different pattern for the off-axis PSF (e.g., measured versus fitted model)
between the data generation process and the training step does not lead to a
significant drop in the detection performance. At this intermediate stage,
each training sample $\widebar{\bm{r}}^{\left[s\right]}$ is obtained by:

$\widebar{\bm{r}}^{\left[s\right]}=\mathbf{P}^{\left[s\right]}\,\bm{r}+\sum\limits_{p=1}^{P^{\left[s\right]}}\alpha_{p}^{\left[s\right]}\,\bm{h}\left(\phi_{p}^{\left[s\right]}\right)\,,$
(9)

where $\mathbf{P}$ is an operator performing the random temporal permutation
of the images of $\bm{r}$ and $\bm{h}\left(\phi_{p}\right)$ in
$\mathbb{R}^{N\times T}$ represents the spatio-temporal contribution of a
synthetic source centered at location $\phi_{p}$ on a reference image at time
$t_{\text{ref}}$, see Sect. 1. The number of sources $P^{\left[s\right]}$,
their contrasts $\\{\alpha_{p}^{\left[s\right]}\\}_{p=1:P^{\left[s\right]}}$
and their initial locations
$\\{\phi_{p}^{\left[s\right]}\\}_{p=1:P^{\left[s\right]}}$ are free
parameters. In practice, the number $P^{\left[s\right]}$ of injected sources
in each training sample is drawn uniformly in $\llbracket 1;10\rrbracket$.
This setting represents a realistic scenario since we expect a few faint
point-like sources in the field of view. We denote by $P$ the total number of
injected sources over the $S$ training sets, i.e.
$P=\sum_{s=1}^{S}P^{\left[s\right]}$. The initial locations
$\\{\phi_{p}\\}_{p=1:P}$ of the injected sources are drawn uniformly per
angular separation (i.e., the number of sources is in average constant per
angular annuli). The range of injected flux is a critical choice. For
instance, if the lower bound is too low, class overlap can occur and the model
is not able to discriminate between sources and the nuisance component leading
to a high level of false alarms. In the opposite case, if the upper bound is
too low, evident bright sources will not be detected since there are no
similar examples in the training set. In practice, we set the injection range
in an unsupervised fashion. We train our model on sources which are
challenging to detect with other methods: the contrast
$\\{\alpha_{p}\\}_{p=1:P}$ of the injected sources is drawn uniformly in
$\left[3\widehat{\sigma}_{\phi_{p}}^{\texttt{PACO}};12\widehat{\sigma}_{\phi_{p}}^{\texttt{PACO}}\right]$
where $\widehat{\sigma}_{\phi_{p}}^{\texttt{PACO}}$ is the 1-sigma contrast
reached by PACO at location $\phi_{p}$. This setting covers both sources that
are detectable above the standard $5\sigma$ detection confidence and sources
whose detection confidence remains below the $5\sigma$ detection limit reached
by PACO. In practice, we found that this setting is suitable to detect both
faint sources and bright sources without generating large number of false
alarms.

As the pre-processing is an expensive procedure and becomes the bottleneck
during online data generation, we adopt a local update strategy to reduce its
computational cost. Prior to the injection of synthetic sources, the whole
dataset is pre-processed, i.e. centered and spatially whitened, see Fig. 1
(b). We denote by $\widetilde{\bm{r}}$ the pre-computed cube. After each batch
$s$ of injections, the set $\mathbb{S}^{[s]}$ of locations impacted by the
signal of the $P^{[s]}$ sources is determined. Outside $\mathbb{S}^{[s]}$, the
pre-processed images are obtained from the temporal permutation of
$\widetilde{\bm{r}}$. Inside $\mathbb{S}^{[s]}$, the statistics of the
nuisance component are updated given the contamination of the $P^{[s]}$
injected sources, and the pre-processed images are updated with these refined
statistics. At this intermediate stage, each training sample
$\widecheck{\bm{r}}^{[s]}$ is obtained by:

$\widecheck{\bm{r}}_{n}^{[s]}=\begin{cases}\mathbf{W}_{n}\left(\widebar{\bm{r}}_{n}^{[s]}\right),&\text{for
}n\in\mathbb{S}^{[s]}\,\cap\,\mathbb{P}\,,\\\
\widetilde{\bm{r}}_{n},&\text{for
}n\in\mathbb{P}-\mathbb{S}^{[s]}\,\cap\,\mathbb{P}\,.\end{cases}$ (10)

This dual strategy, illustrated by Fig. 2, is applied to prevent any detection
bias (i.e., an over-estimation of the actual detection performance of the
proposed algorithm) since we have shown in previous work on the PACO algorithm
(Flasseur et al., 2018b) that the statistics of the nuisance component can
suffer from a (slight) bias, in particular at short angular separations and/or
when the total amount of parallactic rotation is low. This slight potential
bias is due to the contamination of a source whose signal is partly encoded
both in the mean and in the spatial covariances of the nuisance component.

Figure 2: Schematic illustration of the local update of the pre-processing
embedded in the training step of the proposed detection algorithm. (a)
Illustration of the sets $\mathbb{P}$ and $\mathbb{S}^{[s]}$ for a given
training set $s$ with three injected sources displayed in orange. (b)
Illustration of the computation of the pre-processed observations
$\widecheck{\bm{r}}$ in the presence of injected sources from pre-computed
(i.e., before injections) centered and whitened observations
$\widetilde{\bm{r}}$. The spatial scale is not respected for the purpose of
illustration.

Finally, the intermediate images of each training sample are derotated with
the opposite of the parallactic angles so that signal of the synthetic sources
are spatially co-aligned along the temporal axis:

$\mathchoice{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{3.8315pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\displaystyle\bm{r}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{3.8315pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\textstyle\bm{r}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{2.68205pt}{2.72223pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptstyle\bm{r}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{1.91574pt}{1.94444pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptscriptstyle\bm{r}\hss$\crcr}}}\limits}^{[s]}=\mathbf{D}\,\widecheck{\bm{r}}^{[s]}\,,$
(11)

where $\mathbf{D}$ is a derotation operator. This derotation step is mandatory
to perform a semantic segmentation with the CNN we consider (see Sect. 2.2.2)
given the limited spatial extent of its receptive field.

The binary ground-truth segmentation map $\bm{y}^{[s]}$ is obtained by setting
to 1 circular areas centered at the locations $\\{\phi_{p}\\}_{p=1:P^{[s]}}$
of the $P^{[s]}$ injected sources. Other regions of $\bm{y}^{[s]}$ are set to
0. The radius of the circles is set to the full width at half maximum of the
off-axis PSF, which corresponds to the expected spatial extent of an
exoplanetary signature in the data.

A schematic summary of the construction of the training set is given in Fig. 1
(c).

####  2.2.2 Model and architecture

Deep convolutional neural networks reach state-of-art performance in solving
pixel-wise classification tasks for various imaging fields including
microscopy, astronomy, medical imaging or remote sensing. A large variety of
model architecture has been studied in the literature (e.g., auto-encoder
(Badrinarayanan et al., 2017), VGG (Simonyan & Zisserman, 2015), ResNet (He et
al., 2016)) and their performance often rely on an intricate trade-off between
model complexity, the amount of data available for training, and their
fidelity with the data used at test time. A common feature of some classical
deep architectures is to encode the diversity of the training samples in a low
dimensional subspace by transforming the network input with a cascade of
convolution and downsampling operations. Starting from this latent
representation, the initial image size is retrieved by a decoder performing
reverse transformations with a cascade of deconvolution and upsampling
operations.

We also based our model on the above-mentioned category of architectures. We
chose a U-Net (Ronneberger et al., 2015) with a ResNet18 (He et al., 2016) as
encoder backbone ($\simeq$11 millions of free parameters), which is an
architecture widely used for image segmentation. Its residual connections
preserve of the input’s spatial information along the cascade of convolution
and downsampling operations thanks to a direct mapping of the output of each
layer of the compression arm into the corresponding layer of the decompression
arm. We use the architecture implemented in the SMP package555The SMP package
containing the network architecture used in this paper is available at
https://github.com/qubvel/segmentation_models.pytorch.. The encoder and
decoder parts are formed by four blocks, each one is composed by a series of
convolution layers, batch normalization layers, rectified linear unit (ReLU)
activation, and max pooling layers. The final layer of the network has a
sigmoid activation function to produce a detection map
$\widehat{\bm{y}}\in\left[0;1\right]^{M}$. The detailed description of the
architecture, the number of parameters per layer, and the input / output
shapes of each layer can be found at the above mentioned link.

The network weights are trained from scratch with the samples generated with
the procedure described in Sect. 2.2.1. Initial weights are drawn uniformly
through a He-Kaiming distribution (He et al., 2015). The SMP package also
provides pre-trained weights. Pre-training is performed with the ImageNet
dataset (RGB conventional images) either in a supervised, semi-supervised, or
weakly-supervised learning fashion (Yalniz et al., 2019). In case of pre-
training, the weights of the first convolutional layer are replicated in order
to match the $T$-depth of our inputs. We compared all of these strategies
against a supervised learning from scratch with our custom training set (see
Sect. 2.2.1). We found similar performance with all approaches, and opted for
an end-to-end learning.

####  2.2.3 Loss function and accuracy metrics

The design of loss function used for optimizing the network weights at
training time is driven by three criteria: (i) handling with the strong class
imbalance (the number of background pixels being much larger than the number
of pixels from the sources), (ii) being computationally efficient, (iii)
matching the astrophysical goals (i.e., having a measure close to a detection
accuracy score). We compare losses classically used for semantic segmentation,
such as the binary cross-entropy (BCE), $\ell_{1}$ and $\ell_{2}$ norms, mean
square error and Hinge loss. We have also compared losses based on a
similarity measure such as the Dice score (Milletari et al., 2016) and hybrid
losses combining at least two individual loss measurements (e.g., BCE with
Dice score). Our experiments have consistently shown better performance with
Dice-based scores. We selected the Dice2 loss (the 2 means for two classes),
first introduced for biomedical imaging segmentation with very unbalanced
classes (Sudre et al., 2017; Wang et al., 2020). Given a set of ground-truth
and predicted detection maps $\\{\bm{y}^{[s]};\widehat{\bm{y}}^{[s]}\\}$, the
Dice2 score is defined by:

$\displaystyle\mathcal{L}_{\text{detect.}}\left(\bm{y}^{[s]},\widehat{\bm{y}}^{[s]}\right)\,=\,$
$\displaystyle
1-\underbrace{\frac{\sum\limits_{m=1}^{M}(1-\bm{y}_{m}^{[s]})(1-\widehat{\bm{y}}_{m}^{[s]}+\epsilon)}{\sum\limits_{m=1}^{M}2-\bm{y}_{m}^{[s]}-\widehat{\bm{y}}_{m}^{[s]}+\epsilon}}_{\text{background
error}}\,\vspace{1mm}$
$\displaystyle-\,\underbrace{\frac{\sum\limits_{m=1}^{M}\bm{y}_{m}^{[s]}\,\widehat{\bm{y}}_{m}^{[s]}+\epsilon}{\sum\limits_{m=1}^{M}\bm{y}_{m}^{[s]}+\widehat{\bm{y}}_{m}^{[s]}+\epsilon}}_{\text{source
error}}\,,$ (12)

where $\epsilon$ is a minimum-value smoothing and stability parameter added to
avoid division by zero. Targeted loss property (i) is satisfied since errors
in the source and background areas are penalized equally regardless the
relative occurrence of these two classes in $\bm{y}^{[s]}$. Property (ii) is
also satisfied, and we illustrate numerically in the two following paragraphs
that property (iii) is also reached.

Figure 3: Example of training and validation results. (a) Examples of
detection maps obtained at validation time for the best epoch (number 20). (b)
Evolution of the loss function at training and validation time, as well as the
evolution of the F1R accuracy metric at validation time. (c) ROCs representing
the TPR, FDR and F1R as a function of the prescribed detection threshold
$\tau$ for the best epoch (number 20, symbolized by the gray vertical line in
panel (b)). Dataset: HD 95086 (2015-05-05), see Sect. 4 for the recording
logs.

At validation time, we evaluate the ability of the model to detect point-like
sources while simultaneously avoiding false alarms at best as possible. In
other words, we aim to obtain a model obeying a precision-recall trade-off.
For a predicted detection map $\widehat{\bm{y}}^{[s]}$ in
$\left[0;1\right]^{M}$ thresholded at $\tau$ in $\left[0;1\right]$, we count
the number of true positives (TP, i.e., true detections), false positives (FP,
i.e., false alarms) and false negatives (FN, i.e., missed detections).
Following standard practice in direct imaging (see e.g., Flasseur et al.
(2018b); Gonzalez et al. (2018); Cantalloube et al. (2020)), detections are
treated as blobs of one resolution element radius which corresponds to the
expected spatial extent of an exoplanetary signature in the data. From TP, FP,
and FN, we derive the true positive rate (TPR, i.e., the recall), the false
discovery rate (FDR, i.e., the precision), and the F1R score, which is the
harmonic mean between TPR and FDR, that we use as an overall measure of the
precision-recall trade-off:

$\begin{cases}\text{TPR}=\frac{\text{TP}}{\text{TP}+\text{FN}}\in[0;1]\vspace{1mm}\,,\\\
\text{FDR}=\frac{\text{FP}}{\text{FP}+\text{TP}}\in[0;1]\vspace{1mm}\,,\\\
\text{F1R}=\frac{2}{\frac{1}{\text{TPR}}+\frac{1}{\text{FDR}}}=\frac{2\text{TP}}{2\text{TP}+\text{FN}+\text{FP}}\in[0;1]\,.\end{cases}$
(13)

From TPR, FDR, and F1R, receiver operating curves (ROCs; Kay (1993)) are
built. ROCs are obtained by evaluating the figures of merit defined in Eq.
(13) as a function of the detection threshold $\tau$. Finally, the area under
the curve (AUC) for the F1R score is computed as an aggregate score of the
model performance (best when close to 1). Gonzalez et al. (2016, 2018);
Flasseur et al. (2018b); Dahlqvist et al. (2020); Cantalloube et al. (2020);
Daglayan et al. (2022) exemplified the relevance of ROCs in high-contrast
imaging to derive an aggregate measurement of the targeted precision-recall
trade-off.

Figure 3(a) displays some examples of detection maps obtained at validation
time for the best validation epoch. These maps illustrate qualitatively the
ability of our model to detect synthetic sources while simultaneously avoiding
false alarms. Figure 3(b) shows the evolution of the empirical risk (see Eq.
(12)) at training and validation time as well as the evolution of the F1R
accuracy metric (see Eq. (13)) at validation time. The loss function does not
exhibit significant discrepancy between training and validation steps and the
convergence is reached in a few epochs. Besides, the accuracy score is high
and well anti-correlated with the loss. This latter observation illustrates
that the loss function is a satisfactory estimate of the overall accuracy
metric (see targeted property (iii) at the beginning of Sect. 2.2.3). Finally,
Fig. 3(c) gives an illustration of validation ROC obtained for the best epoch
(symbolized by a gray vertical line in Fig. 3(b)). The validation ROC confirms
the good ability of the trained model to discriminate (synthetic) point-like
sources from the nuisance component.

####  2.2.4 Implementation details

Pairs of samples $\\{\mathchoice{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{3.8315pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\displaystyle\bm{r}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{3.8315pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\textstyle\bm{r}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{2.68205pt}{2.72223pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptstyle\bm{r}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{1.91574pt}{1.94444pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptscriptstyle\bm{r}\hss$\crcr}}}\limits}^{[s]};\bm{y}^{[s]}\\}$
are generated on the fly at training and evaluation time following the
procedure described in Sect. 2.2.1. To avoid over-fitting, each realization
$s$ is unique with no repetition for the different epochs. The notion of epoch
is used only as a way to evaluate regularly the performance of the model with
the validation procedure described in Sect. 2.2.3, and also to manage the
learning rate scheduler. The optimization of the network weights is performed
with an iterative stochastic gradient-descent strategy on mini-batches of
samples. It means that, at each iteration, the model weights are updated in
the opposite direction to the gradient of the loss. The loss is evaluated from
the current mini-batch of samples with respect to the model weights. Since we
work on series of $T$ images, with $T$ typically lying between 50 and 100
images, the batch-size is fixed at 1 given memory constraints. This setting
does not affect the overall performance of the method and only requires to
perform more iterations to reach convergence since the cost function is quite
noisy, see Fig. 3(a). Even under this setting, the convergence is typically
reached in a few epochs, see Sect. 2.2.3 and Fig. 3(b). In practice, for each
training epoch, $S=100$ pairs of samples
$\\{\mathchoice{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{3.8315pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\displaystyle\bm{r}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{3.8315pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\textstyle\bm{r}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{2.68205pt}{2.72223pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptstyle\bm{r}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{1.91574pt}{1.94444pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptscriptstyle\bm{r}\hss$\crcr}}}\limits}^{[s]};\bm{y}^{[s]}\\}$
are generated and fed sequentially as input of the network. For each
validation epoch, $S$ is fixed at 10 given the computational burden required
to build ROCs representing the F1R score as a function of the detection
threshold $\tau$. The training process stops when the accuracy metric
$\text{AUC}_{\text{F1R}}$ (i.e., AUC under ROC representing the F1R score as a
function of the detection threshold $\tau$ at validation time) evolves in less
than 2% during the ten previous epochs. The model optimization is performed
with the adaptive gradient algorithm AMSGrad (Reddi et al., 2019) which is a
variant of the Adam (Kingma & Ba, 2014) optimizer with a longer-term memory of
past gradients. The parameters of the optimizer and of the scheduler have been
fine tuned on two datasets and are kept constant for all our experiments. In
practice, we observed that the optimized values are quite robust with respect
to the dataset diversity. The weight decay is fixed at $10^{-5}$ and the
initial learning rate is set to $10^{-3}$ with a regular decrease by 10% every
10 epochs. The optimization of the network weights is performed with the high-
performance deep learning library PyTorch (Paszke et al., 2019) on GPUs server
with NVIDIA system equipped with either Tesla V100 or GTX 1080 Ti cards. The
pre-processing step is highly parallelized and has a double implementation so
that it can performed either on CPUs or on GPUs depending on the number of
available CPUs cores and on the server specifications.

##  3 Source characterization algorithm

Once sources have been detected, they can be characterized in terms of
astrometry and of photometry. In this section, we present a new method based
on supervised deep learning to estimate the photometry of detected point-like
sources. The sub-pixel estimation of the astrometry is not addressed in this
paper because it requires a sub-pixel estimation of the detection criterion as
well as statistical guarantees on its significance. These specific
developments are left for future work, and the proposed characterization
algorithm can be use to estimate the photometry at the pixel barycenter of
(candidate) sources revealed by the detection stage presented in Sect. 2. As
in Sect. 2, we successively discuss the pre-processing stage, the
formalization of the problem as a regression task, the model and the
underlying architecture, the metric we use for training and evaluation, and
some implementation details.

###  3.1 Pre-processing aspects

Figure 4: Influence of the whitening of the spatial correlations during the
pre-processing step of the characterization algorithm. ARE (see Sect. 3.2.3
for the definition of the metric) on the estimated photometry is reported as a
function of the angular separation, with and without whitening of the spatial
correlations. The performance of PACO are displayed as a purpose of
comparison. The results are averaged azimuthally for 40,000 sources of flux
drawn uniformly between $1\times 10^{-6}$ and $3\times 10^{-5}$. The known
real sources were masked out and were not considered. Dataset: HD 95086
(2015-05-05), see Sect. 4 for the recording logs.

We adopt a simple patch-based approach, in which we predict the flux of a
putative source from a unique (reduced) patch centered on the approximate
location of the source. We propose to parameterize the mapping between the
input patch and the flux with a CNN trained in a supervised fashion, from the
dataset of interest (i.e., the underlying model is data-dependent, as for the
detection stage of this paper).

The dataset is first reduced to a single frame, from which the input patches
are extracted. This pre-processing consists in three steps. First, the
temporal mean is computed and subtracted pixel-wise in order to remove most of
the quasi-static speckles. Second, the dataset is derotated by the opposite of
the parallactic angles to co-align the signal of the sources. Third, the
dataset is averaged temporally, resulting in a single averaged frame. This
last step allows to obtain an efficient training procedure as it reduces the
size of the input data by a factor $T$. Besides, we observed empirically that
this step is beneficial to improve the estimation accuracy as it acts as a
simple denoiser: the source signal is constant along time, while residual
speckles are not quasi-static after cube derotation, thus canceling out.
Keeping the notation introduced in Sect. 2, these operations transform a given
(intermediate) training dataset $\widebar{\bm{r}}$ in $\mathbb{R}^{N\times T}$
with injected synthetic sources as follows:

$\begin{cases}\widecheck{\bm{r}}_{n}=\widebar{\bm{r}}_{n}-\widehat{\bm{m}}_{n}\,,\forall
n\in\mathbb{P}\,,\text{(step 1)}\,,\\\
\mathchoice{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{3.8315pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\displaystyle\bm{r}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{3.8315pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\textstyle\bm{r}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{2.68205pt}{2.72223pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptstyle\bm{r}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{1.91574pt}{1.94444pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptscriptstyle\bm{r}\hss$\crcr}}}\limits}=\frac{1}{T}\sum\limits_{t=1}^{T}\left[\mathbf{D}\,\widecheck{\bm{r}}\right]_{t}\,,\text{(steps
2 and 3)}\,.\end{cases}$ (14)

In this framework, and unlike the detection stage, we do not apply a whitening
of the spatial correlations at step 1. Indeed, we observed empirically that
whitening the dataset between steps 1 and 2-3 degrades the performance of our
model, as illustrated by Fig. 4. This is expected as the whitening distorts
the shape and the norm of the source signal, thus hampering the recovery of
its flux. The detection algorithm is not subject to this constraint as its
task is to determine whether a source is present or not, regardless of its
flux. After pre-processing, squared patches
$\mathchoice{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{4.02501pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\displaystyle\bm{p}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{4.02501pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\textstyle\bm{p}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{2.8175pt}{2.72223pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptstyle\bm{p}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{2.0125pt}{1.94444pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptscriptstyle\bm{p}\hss$\crcr}}}\limits}^{[p]}\in\mathbb{R}^{J}$
are finally extracted around the location of each injected synthetic source
$p$ during the training and validation steps, or around each (candidate) real
point-like source $p$ at inference time. The patch size $J$ is an hyper-
parameter whose setting is discussed in more details in Sect. 3.2.4.

Figure 5: Influence of the patch size $J$ in the characterization algorithm.
ARE (see Sect. 3.2.3 for the definition of the metric) on the estimated
photometry is reported as a function of the angular separation, for patches
between $J=15^{2}$ and $J=31^{2}$ pixels area. The results are averaged
azimuthally for 40,000 sources of flux drawn uniformly between $1\times
10^{-6}$ and $3\times 10^{-5}$. The known real sources were masked out and
were not considered. Dataset: HD 95086 (2015-05-05), see Sect. 4 for the
recording logs.

###  3.2 Regression by supervised deep learning

####  3.2.1 Construction of the training basis

As with the detection procedure of Sect. 5, we resort to massive injections of
synthetic sources with various fluxes to build our training basis.

Prior to the injections, the first step consists in masking any real and/or
synthetic detected sources that we aim to estimate the photometry at inference
time. This operation prevents, as best as possible, data leakages between
training and test sets so that training patches do not contain any pixel from
patches considered at inference time. In practice, source masking is performed
by replacing, for each temporal frame, the local area impacted by the signal
of the sources of interest by their pixel-wise temporal mean. Like for the
generation of the training basis of the detection stage, we also apply, as a
data-augmentation strategy, a random permutation of the temporal frames prior
to the construction of a training set $s$ from which training samples with
injected sources are extracted.

We then build each training set $s$ by injecting a dozen of synthetic sources,
with a relative flux (i.e., exoplanet to star contrast) ranging from $1\times
10^{-6}$ to $3\times 10^{-5}$ with respect to the flux of the host star. Given
the current instrumental and processing performance in direct imaging, this
setting corresponds to sources with relatively low flux, for which the
estimation is usually the most flawed. This operation range can be easily
modified in the algorithm to characterize (expected) sources with a lower or
higher contrast level, if needed. The contrast of synthetic sources is drawn
uniformly in the above-mentioned range, regardless of the angular separation.
After injection, each training set is pre-processed and the input patches are
extracted following the method described in Sect. 3.1. This procedure is
repeated to get a prescribed number $P$ of training patches
$\\{\widecheck{\bm{p}}^{[p]}\in\mathbb{R}^{J}\\}_{p=1:P}$. The number of
training patches, hence the total number $P$ of injected synthetic sources, is
an additional hyper-parameter whose setting is discussed in more details in
Sect. 3.2.4.

####  3.2.2 Model and architecture

Table 2: Architecture of the proposed CNN for source characterization. The shapes of the layers are indicated for a unit batch size. layer | shape
---|---
Input | $1\times 31\times 31$
Conv2D + ReLU | $128\times 25\times 25$
Conv2D + ReLU | $128\times 21\times 21$
Conv2D + ReLU | $256\times 17\times 17$
MaxPooling | $256\times 1\times 1$
DenseLayer + ReLU | $256\times 1\times 1$
DenseLayer | $1\times 1\times 1$

We built a custom network based on VGG, an architecture initially proposed for
image classification (Simonyan & Zisserman, 2015). The underlying model has
1.2 million of free parameters, and its detailed architecture is described in
Table 2. We use a stride of one for all convolutional layers. We have also
tested alternative models, both with a deeper and shallower architecture, all
leading to worst estimation performance than the selected one, especially at
short angular separations.

From an input patch $\mathchoice{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{4.02501pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\displaystyle\bm{p}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{4.02501pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\textstyle\bm{p}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{2.8175pt}{2.72223pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptstyle\bm{p}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{2.0125pt}{1.94444pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptscriptstyle\bm{p}\hss$\crcr}}}\limits}^{[p]}\in\mathbb{R}^{J}$,
the network produces a single scalar
$\widehat{\alpha}^{[p]}\in\mathbb{R}_{+}$, representing the estimated source’s
photometry.

####  3.2.3 Loss function and accuracy metric

Our choice of the loss function $\mathcal{L}_{\text{carac.}}$ is driven by the
two following criteria: (i) being computationally efficient, and (ii) matching
the astrophysical goals. We chose the absolute relative error (ARE) between
the ground-truth and the predicted photometry, which is a classical loss for
regression problems:

$\mathcal{L}_{\text{carac.}}\left(\alpha^{[p]},\widehat{\alpha}^{[p]}\right)\,\,\,(\%)=100\times\frac{\big{|}\alpha^{[p]}-\widehat{\alpha}^{[p]}\big{|}}{\alpha^{[p]}}\,.$
(15)

The ARE has the advantage of giving the same contribution to each individual
source $p$ regardless of its flux, when it is averaged over multiples ones.
Since this metric is computationally very efficient, we also use it to
evaluate the overall performance of the method at validation time.

####  3.2.4 Implementation details

In this section, we successively discuss the setting of the patch size, the
number of training sources, the sampling strategy of injected synthetic
sources, and some optimization aspects.

Due to the pre-processing described in Sect. 3.1, part of the signal of the
sources can be encoded in the temporal mean that is subtracted to the full
frames in order to attenuate the quasi-static speckles. As a result, a
(negative) contribution, taking the form of an arc, can spread out along the
trajectory of the sources in the reduced frame. This well-known phenomenon in
direct imaging is usually referred as self-subtraction. As a consequence, we
can expect that the performance of the predictor would increase with the patch
size $J$, to be able to capture the (extended) signature of sources induced by
self-subtraction. Besides, increasing the patch size increases the context
(i.e., local realizations of the nuisance component) that could be beneficial
to unmix the different contributions. As shown by Fig. 5, increasing the patch
size is indeed beneficial as it reduced the mean relative error of estimation.
The gain is more pronounced as the angular separation increases, which could
be due to the larger extent of the self-subtraction signature (as the apparent
motion of sources induced by ADI increases with the angular separation).
However, large patches are not convenient in the case of adjacent sources, as
both signals will be contained in both input patches. As a trade-off, we chose
a patch size of $J=31^{2}$ pixels, as it encompasses the core of the signal of
the source without being impractical when multiple sources are present in the
field of view.

The number $P$ of synthetic sources (which also corresponds to the number of
patches) used at training time is an additional hyper-parameter obeying a
trade-off. On the one hand, it should be large enough to be representative of
the variety of real sources in terms of flux and locations. On the other hand,
it should be sufficiently small to avoid over-fitting on the training set.
These expected behaviors are confirmed by numerical experiments presented in
Fig. 6. The error at small angular separations increases with $P$ since this
is the area of the field of view the more subject to data leakages when
generating multiple samples from a few tens of pixels only. The overall
performance also degrades when $P$ is too small. Based on this study, we
include $P=40,000$ patches in our training set since it leads to the smallest
ARE averaged over the whole of view.

Concerning the source’s sampling strategy, we are interested at evaluation
time in assessing the performance of the model per angular separation. As
such, it is natural to sample test sources uniformly in the polar coordinates
system. However, polar sampling is detrimental during the training phase, as
pixels at short angular separations would be over-represented in the training
set, leading to an over-fitting of the model in this area. We experimentally
observed that sampling training sources uniformly in the Cartesian coordinates
system reduces significantly over-fitting at short angular separations,
without degrading performance in the rest of the field of view. As a result,
we opt during training for a uniform sampling of the coordinates of synthetic
sources in the Cartesian system.

Figure 6: Influence of the number of training sources $P$ in the
characterization algorithm. ARE (see Sect. 3.2.3 for the definition of the
metric) on the estimated photometry is reported as a function of the angular
separation, for $P\in\\{$400,000; 40,000; 20,000; 4,000$\\}$ sources. The
results are averaged azimuthally for sources of flux drawn uniformly between
$1\times 10^{-6}$ and $3\times 10^{-5}$. The mean ARE (denoted by
$\overline{\text{ARE}}$) averaged over the angular separation is also reported
as an overall measure of the performance. The known real sources were masked
out and were not considered. Datasets: the eleven SPHERE-IRDIS datasets
considered in this work, see Sect. 4 for the recording logs.

Concerning the optimization process, pairs of samples
$\\{\mathchoice{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{4.02501pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\displaystyle\bm{p}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{4.02501pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\textstyle\bm{p}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{2.8175pt}{2.72223pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptstyle\bm{p}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{2.0125pt}{1.94444pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptscriptstyle\bm{p}\hss$\crcr}}}\limits}^{[p]};\alpha^{[p]}\\}$
are pre-generated before training and evaluation. As for the detection stage,
each realization $p$ is unique to limit over-fitting as best as possible.
Unlike the detection stage, the characterization stage operates on a single
patch (instead of the $T$ temporal images), thus allowing to chose a batch
size higher than one to improve the stability and to reduce the computation
time of the optimization process. In practice, the batch size is fixed at
1024. The number of epochs is fixed at 300, which shown to be sufficient to
reach convergence of the network weights for all the considered observations.
For each training epoch, the full set of training samples is fed as input of
the network in a random order. The model optimization is performed with Adam
(Kingma & Ba, 2014) with a learning rate of $10^{-3}$. As for the detection
stage, the optimization of the network weights is done in PyTorch (Paszke et
al., 2019) on either Tesla V100 or GTX 1080 Ti cards. The pre-processing step
is highly parallelized, and it is run on CPUs.

##  4 Results

###  4.1 Datasets description and reduction strategies

Table 3: Observing parameters of ADI sequences from the VLT/SPHERE-IRDIS instrument considered in this paper. Columns are: target name, ESO survey ID, observation date, spectral filter $\lambda$, number $T$ of available temporal frames, total apparent rotation $\Delta_{\text{par}}$ of the field of view, number NDIT of sub-integration exposures, individual exposure time DIT, average coherence time $\tau_{0}$, average seeing, and the first paper reporting analysis of the same data. All the observations are performed with the apodized Lyot coronagraph (Carbillet et al., 2011) of the VLT/SPHERE instrument. ${}^{\text{(a)}}$Since the EIDC aimed to perform a blind benchmark, information that would allow to identify the datasets are not known. ${}^{\text{(b)}}$This dataset was recorded with the star-hopping technique recently available for the VLT/SPHERE instrument (Wahhaj et al., 2021) and its analysis is not reported yet. We do not exploit the dataset associated to the observation of the reference star. The dataset associated to the target star (HD 95086) is processed like all other datasets considered in this paper. ${}^{\text{(c)}}$This dataset is only used for additional experiments conducted in Appendix A.2. Target | ESO ID | Obs. date |  $\lambda$ ($\mathrm{\SIUnitSymbolMicro}\mathrm{m}$) | $T$ |  $\Delta_{\text{par}}$ (°) | NDIT | DIT (s) |  $\tau_{0}$ (ms) | Seeing (”) | Related paper
---|---|---|---|---|---|---|---|---|---|---
VLT/SPHERE-IRDIS observations from the EIDC challenge
IRDIS 1${}^{\text{(a)}}$ | $-^{\text{(a)}}$ | $-^{\text{(a)}}$ | 1.625 | 252 | 40.3 | $-^{\text{(a)}}$ | $-^{\text{(a)}}$ | $-^{\text{(a)}}$ | $-^{\text{(a)}}$ | Cantalloube et al. (2020)
IRDIS 2${}^{\text{(a)}}$ | $-^{\text{(a)}}$ | $-^{\text{(a)}}$ | 1.593 | 80 | 31.5 | $-^{\text{(a)}}$ | $-^{\text{(a)}}$ | $-^{\text{(a)}}$ | $-^{\text{(a)}}$ | Cantalloube et al. (2020)
IRDIS 3${}^{\text{(a)}}$ | $-^{\text{(a)}}$ | $-^{\text{(a)}}$ | 1.593 | 228 | 80.5 | $-^{\text{(a)}}$ | $-^{\text{(a)}}$ | $-^{\text{(a)}}$ | $-^{\text{(a)}}$ | Cantalloube et al. (2020)
VLT/SPHERE-IRDIS observations
HD 95086 | 095.C-0298(A) | 2015-05-05 | 2.110 | 52 | 18.2 | 4 | 64 | 2.3 | 0.89 | Chauvin et al. (2018)
HD 95086 | 1100.C-0481(E) | 2018-01-05 | 2.110 | 70 | 41.0 | 10 | 96 | 7.8 | 0.32 | Desgrange et al. (2022)
HD 95086 | 106.21VL.001 | 2021-03-11 | 2.110 | 104 | 41.4 | 2 | 32 | 7.0 | 0.77 | $-^{\text{(b)}}$
HIP 88399 | 095.C-0298(A) | 2015-05-10 | 1.593 | 46 | 34.3 | 4 | 64 | 1.2 | 1.05 | Langlois et al. (2021)
HIP 88399 | 097.C-0865(A) | 2016-04-16 | 1.593 | 54 | 37.3 | 5 | 64 | 2.0 | 1.45 | Langlois et al. (2021)
HIP 88399 | 1100.C-0481(F) | 2018-04-11 | 1.593 | 40 | 31.9 | 10 | 96 | 5.5 | 0.74 | Langlois et al. (2021)
HD 131399 | 095.C-0389(A) | 2015-06-12 | 2.110 | 92 | 36.7 | 6 | 16 | 1.9 | 0.90 | Wagner et al. (2016)
HD 131399 | 296.C-5036(A) | 2016-05-07 | 2.110 | 56 | 39.5 | 7 | 32 | 3.6 | 0.98 | Wagner et al. (2016)
HIP 65426 | 198.C-0209(E) | 2017-02-09 | 2.110 | 55 | 49.1 | 4 | 64 | 4.4 | 0.82 | Chauvin et al. (2017)
HIP 65426 | 1100.C-0481(G) | 2018-05-13 | 2.110 | 40 | 31.7 | 10 | 96 | 4.3 | 0.81 | Cheetham et al. (2019)
HIP 72192 | 095.C-0389(A) | 2015-06-11 | 2.110 | 96 | 17.3 | 6 | 16 | 1.9 | 1.03 | Flasseur et al. (2018b)
HR 8799${}^{\text{(c)}}$ | 095.C-0298(C) | 2015-07-04 | 2.110 | 112 | 17.9 | 8 | 32 | 2.3 | 0.94 | Langlois et al. (2021)

For our comparative analysis, we have selected 15 datasets recorded with the
SPHERE-IRDIS instrument.

Three of the 15 SPHERE-IRDIS datasets were extracted from the exoplanet
imaging data challenge (EIDC) initially designed to ground the detection
performance of existing post-processing algorithms for high-contrast imaging
(Cantalloube et al., 2020). These datasets are used as a sanity check to
assess the ability of the proposed algorithm to detect injected sources at
moderate levels of contrast.

To study in more details the performance of the proposed method, we selected
12 additional datasets, mostly part from the SHINE survey of the SPHERE-IRDIS
instrument (Desidera et al., 2021; Langlois et al., 2021; Vigan et al., 2021).
They were obtained by the observation of the following stars:
– HD 95086, a A8III type star of the Carina constellation, hosting an
exoplanet discovered by direct imaging with the SPHERE instrument (Rameau et
al., 2013a, b). Ten known background point sources are also within the SPHERE-
IRDIS field of view. In addition, based on the analysis of PACO and deep PACO
detection maps, we have identified two additional (unbounded) candidate point-
like sources. Given their unknown status, we exclude these sources from our
general analysis (i.e., there are not considered as true detections or false
alarms). We briefly discuss there status in Sect. 4.2.3.
– HIP 88399, a F6V type star of the Vela constellation, without known bounded
exoplanet. However, six faint background sources fall within the SPHERE-IRDIS
field of view.
– HD 131399 A, a A1V type star of a triple system located in the Centaurus
constellation, with a known faint point source (HD 131399 Ab) discovered by
direct imaging (Wagner et al., 2016). While first supposed to be an exoplanet,
follow-up analysis of its astro-photometry show that HD 131399 Ab is more
likely a background brown dwarf (Nielsen et al., 2017). Besides, a bright
background star falls within the SPHERE-IRDIS field of view.
– HIP 65426, a A8III type star of the Carina constellation, hosting an
exoplanet discovered by direct imaging with the SPHERE instrument (Chauvin et
al., 2017).
– HIP 72192, a A0V type star of the Lupus constellation, without known bounded
exoplanet. However, two faint background sources fall within the SPHERE-IRDIS
field of view.
Considered SPHERE-IRDIS observations were obtained in H2-H3 (i.e.,
$\lambda_{\textit{H2}}=1.593\,\mathrm{\SIUnitSymbolMicro}\mathrm{m}$,
$\lambda_{\textit{H3}}=1.667\,\mathrm{\SIUnitSymbolMicro}\mathrm{m}$) or in
K1-K2 (i.e.,
$\lambda_{\textit{K1}}=2.110\,\mathrm{\SIUnitSymbolMicro}\mathrm{m}$,
$\lambda_{\textit{K2}}=2.251\,\mathrm{\SIUnitSymbolMicro}\mathrm{m}$) dual
spectral bands. The main parameters of each observation are summarized in
Table 3. The diversity in the experienced observing conditions is quite
representative of the SPHERE observations.

The raw observations were pre-reduced666In this paper, we use the term pre-
reduction to stand for the extraction, mapping, and correction of the raw
data. We use the term pre-processing to stand for the centering and whitening
of the pre-reduced observations with the statistical model detailed in Sect.
2.1, that serve as inputs of the supervised deep learning stage detailed in
Sect. 2.2. with the DRH pipeline (Pavlov et al., 2008) of the SPHERE
instrument, which performs thermal background subtraction, flat-field
correction, anamorphism correction, compensation for spectral transmission,
flux normalization, bad pixels identification and interpolation, frame
centering, true-North alignment, wavelength calibration, astrometric
calibration, and frame selection. These operations are complemented by custom
routines implemented in the SPHERE data center (Delorme et al., 2017), in
particular to improve bad pixels correction. Finally, the SPHERE data center
combines the pre-reduced observations and delivers the calibrated ADI datasets
we consider in this work. ADI reduction is performed by considering the first
spectral channel (i.e., either at
$\lambda_{\textit{H2}}=1.593\,\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ or at
$\lambda_{\textit{K1}}=2.110\,\mathrm{\SIUnitSymbolMicro}\mathrm{m}$) given
that the contrast is significantly more favorable in this channel than in the
second one.

To ground the performance of the proposed algorithm, we resort to massive
injections of synthetic sources, as done to train the deep model of the
proposed algorithm, see Sect. 2.2.1. Given the simplicity of the signature of
the sought objects (i.e., taking the form of a blob spatially correlated over
only a few pixel width), we did not find significant bias in using the same
injection procedure for the training and evaluation steps. Injections of
synthetic sources in the EIDC benchmark were performed by there authors with
the VIP pipeline.

The detection performance of the proposed method are compared in Sect. 4.2
with the cADI, PCA, and PACO algorithms, (see Sect. 1 for their respective
principle). For cADI, we have re-implemented the original method (Marois et
al., 2006) based on a full-frame estimation of the off-axis PSF and of the S/N
map, i.e., without angular-specific processing. We have also used the refined
implementation of cADI available in the VIP package (Gonzalez et al., 2017),
which includes a protection angle strategy accounting for a minimal field
rotation between successive images when building the off-axis PSF in order to
limit the self-subtraction effect. After computation of the off-axis PSF, a
S/N map is derived by accounting for an annular-based estimation of the noise
in the residual images. We also applied the VIP implementation of the PCA-
based algorithm combined with the same protection angle strategy and the
annular-based computation of the S/N. For PCA reductions, the number of modes
has been optimized in annuli by maximizing the S/N of synthetic sources with
similar ranges of contrast than the ones we consider for our comparisons. The
other parameters of the VIP implementation of cADI and PCA are less critical
and are fixed at pre-set values (Gonzalez et al., 2017). For PACO, we
performed the data reduction with our fully unsupervised processing pipeline
(Flasseur et al., 2018a).

Concerning the photometry estimation, we compare in Sect. 4.3 the performance
of the proposed method with PACO and the VIP implementation of the PCA. PACO
parameters are estimated automatically in a data-driven fashion. In a
nutshell, the characterization of a point-like source $p$ is performed by a
joint estimation of (i) the statistics (i.e., mean $\bm{m}_{n}$ and
covariances $\mathbf{C}_{n}$) of the nuisance component $\bm{f}$, and (ii) of
the photometry (i.e., flux $\alpha_{p}$) of the given source, see Eq. 1. For
the PCA, we resort to a similar procedure than for the detection step to set
the parameters. In particular, the number of modes is optimized for each
injected source to be characterized by maximizing its S/N. Once the setting
fixed, the flux of a given source is estimated by minimizing the residuals
through the injection of negative fake companions (Wertz et al., 2017). This
is performed with a two steps procedure, as recommended in the VIP package:
(i) a first guess estimate is obtained by performing a grid-search, (ii) a
local optimization is performed with a Nelder-Mead simplex algorithm (Nocedal
& Wright, 1999). Given computational constraints (largely dominated by the
PCA), the exact (known) sub-pixel location $\phi_{p}$ of each injected source
$p$ is provided to the different algorithms (i.e., it is not optimized), and
the photometry is estimated at this ground-truth position. When estimating the
photometry of real sources, both the astrometry and the photometry are
optimized by the different algorithms.

###  4.2 Detection results

####  4.2.1 Detection of known real sources

Figure 7: Detection maps obtained with the selected algorithms (see Sect.
4.1). Sources are classified as true, missed and false detections. The two
additional candidate point-like sources whose identification is discussed in
Sect. 4.2.3 are not classified. The detection threshold is set to $\tau=5$ for
cADI, cADI (VIP), PCA (VIP) and PACO. It is set to $\tau=0.5$ for the proposed
algorithm. The light blue line represents the sensor field of view while the
dashed blue line represents the extended field of view on which the detection
can be performed due to the apparent rotation of the field induced by ADI.
Dataset: HD 95086 (2015-05-05), see Sect. 4 for the observation logs.

Figure 8: Same caption than Fig. 8. Zoom near the host star. Dataset: HD 95086 (2015-05-05), see Sect. 4 for the observation logs. Figure 9: ROCs showing the TPR as a function of the FDR for real sources. Dataset: HD 95086 (2015-05-05), see Sect. 4 for the observation logs. Table 4: Mean results of AUC for ROCs giving the TPR as a function of the FDR. The scores are averaged over the eleven SPHERE-IRDIS datasets considered in this paper, see Sect. 4 for the observation logs. Only the 59 known real sources present in these datasets were considered. Figure 7 of (Flasseur et al., 2022) display the corresponding ROC from which these mean results were aggregated. sep. (”) | cADI | cADI (VIP) | PCA (VIP) | PACO | proposed
---|---|---|---|---|---
$\left[0;9\right]$ | 0.38 | 0.72 | 0.66 | 0.88 | 0.95

A first classical sanity test to evaluate the detection performance of a post-
processing algorithm is to study qualitatively its ability to re-detect real
known sources initially detected with different algorithms, and possibly from
different datasets. We present detailed results for one dataset (HD 95086,
2015-05-05, see Table 3) selected among the eleven SPHERE-IRDIS observations
we consider because HD 95086 is the star having the larger number of known
real sources in the SPHERE-IRDIS field of view. Results obtained for the ten
other datasets are given in supplementary material. Figures 8 and 8 give
detection maps produced with the five tested algorithms. The detection
threshold is set to $\tau=5$ for the algorithms producing a S/N map (i.e.,
cADI, cADI (VIP), PCA (VIP), PACO), and to $\tau=0.5$ for the proposed method
producing a pseudo-probability map. Due to the binary pixel-wise
classification task we consider for the training step of the proposed method
(see Sect. 2.2), its detection map is almost binary (i.e., each pixel value is
close either to 0 or 1) so that the setting of the threshold $\tau$ is quite
flexible. Based on the analysis of the detection maps, PACO and the proposed
method lead to the best qualitative results since there are the only
algorithms able to detect all real known sources without any false alarm in
most of the field of view. With the proposed method, only two false alarms
occur very near the borders of the field of view due to the limited number
(much lower than $T$) of temporal samples available in this area to build a
consistent model of the nuisance. This claim is also supported by the PACO
detection map that displays a few false alarms in the same area.

Figure 9 gives a more quantitative analysis of the previous results through
ROCs representing the TPR as a function of the FDR (see Sect. 2.2.3 and Eq.
(13)) for the same dataset of HD 95086 (2015-05-05). This type of
representation gives a comparison of the precision-recall trade-off reached by
each method, regardless the detection quantity (S/N or pseudo-probability)
they produce. These curves are obtained by counting the number of true
positives (TPs), and false alarms (FAs) for the full range of possible
detection thresholds, i.e. $\tau\in\left[0;1\right]$ for the proposed method,
and
$\tau\in\left[\text{min}(\widehat{\bm{y}});\text{max}(\widehat{\bm{y}})\right]$
for cADI, PCA, and PACO. Table 4 presents averaged results over the eleven
SPHERE-IRDIS datasets we consider in this study, and the detailed scores for
each dataset are given in Table 1 of the supplementary material. These results
illustrate the benefits of the proposed method in terms of precision-recall
trade-off: the AUC under ROC is improved by at least 7$\%$ with respect to the
comparative algorithms.

####  4.2.2 Detection of synthetic sources

Figure 10: Detection maps obtained with the proposed algorithm on the three SPHERE-IRDIS datasets from the EIDC data challenge, see Sect. 4 for the observation logs. Sources are classified as true, missed and false detections. The detection threshold is set to $\tau=0.5$. Table 5: Detection scores: F1R at detection threshold $\tau=0.5$ (the higher, the better), AUC under the ROC representing the TPR as a function of $\tau$ (the higher, the better), and AUC under the ROC representing the FDR as a function of $\tau$ (the lower, the better). Scores reported for PACO are extracted from the EIDC data challenge (Cantalloube et al., 2020), and deep PACO scores are computed with a similar procedure for the three SPHERE-IRDIS datasets from the EIDC. –${}^{\text{(a)}}$Metrics can not be computed since there is no injected source in this dataset. | IRDIS 1 | IRDIS 2 | IRDIS 3 | mean | rank | IRDIS 1 | IRDIS 2 | IRDIS 3 | mean | rank
---|---|---|---|---|---|---|---|---|---|---
|  PACO | deep PACO
F1R | 1.00 | –${}^{\text{(a)}}$ | 1.00 | 1.00 | 1st/22 (on par) | 1.00 | –${}^{\text{(a)}}$ | 1.00 | 1.00 | 1st/23 (on par)
$\textbf{AUC}_{\textbf{TPR}}$ | 1.00 | –${}^{\text{(a)}}$ | 0.93 | 0.97 | 1st/22 (on par) | 1.00 | –${}^{\text{(a)}}$ | 0.96 | 0.98 | 1st/23
$\textbf{AUC}_{\textbf{FDR}}$ | 0.39 | –${}^{\text{(a)}}$ | 0.32 | 0.36 | 6/22 | 0.01 | –${}^{\text{(a)}}$ | 0.01 | 0.01 | 1st/23

In this section, we evaluate quantitatively the detection performance of the
proposed method with ROCs and contrast curves built from massively injected
synthetic sources.

As a first step, we apply the proposed detection algorithm on the three
SPHERE-IRDIS datasets from the public EIDC data challenge (Cantalloube et al.,
2020). The resulting detection maps are given in Fig. 10. Using the same
procedure than Cantalloube et al. (2020) to benchmark 22 post-processing
algorithms (including cADI, PCA, and PACO), we compute the F1R score at the
set detection threshold ($\tau=0.5$), and ROCs for the TPR and the FDR scores,
as defined in Eqs. (13), by varying the detection threshold. The AUC is then
computed from ROCs. Table 5 summarizes the results we obtained with the
proposed algorithm. As a purpose of comparison, we also report the PACO
results that have been published in Cantalloube et al. (2020). It emphasizes
the interesting precision-recall of the proposed approach that performs, on
these datasets, on par with or better than the 22 post-processing algorithms
considered in the EIDC data challenge. However, these results should be taken
with caution since they are based on a few datasets, with only six injected
sources at relatively bright levels of contrast. Besides, several algorithms
(including PACO) are also able to detect the injected sources without any
false alarm in the field of view for a sufficiently large detection threshold.
At this stage, the better performance of deep PACO in terms of false alarms
rejection (quantified by the $\text{AUC}_{\text{FDR}}$ metric) are mostly
explained by the fact that the proposed algorithm produces detection maps with
almost binary values. It can be noted that this effect is also encountered for
all the algorithms of the EIDC data challenge producing the same type of
outputs like RSM (Dahlqvist et al., 2020) or SODINN (Gonzalez et al., 2018).

To ground in details the precision-recall trade-off of the proposed detection
algorithm, it is necessary to build consistent ROCs and contrast curves by
resorting to massive injections of synthetic sources at various levels of
contrast. For that purpose, the most realistic procedure, hereafter called
reference procedure, consists in (i) splitting the whole set of synthetic
sources in small subsets, (ii) injecting synthetic sources of one subset in
the dataset of interest so that injected fake sources mimic the behaviour of
real (possibly unknown) sources, (iii) training the detection model with the
procedure described in Sect. 2, and (iv) applying the trained model to the
dataset containing the synthetic sources injected in step (ii). Steps (ii) to
(iv) are repeated for all subsets of synthetic sources. This procedure
simulates the real situation when we face a new dataset with real unknown
exoplanets that we aim to detect at inference time. Due to the computational
burden of step (iii), repeating this full procedure for all subsets of
synthetic sources is not realistic. To circumvent this issue, we resort to a
proxy procedure: instead of training a different model for each subset of
injected sources, we train a unique model without synthetic sources mimicking
the behaviour of real (possibly unknown) sources. Synthetic sources are
injected a posteriori, i.e., after completing the training step (iii), and the
trained model is applied on the resulting dataset. This proxy procedure leads
to an improvement in terms of algorithmic complexity by a factor equal to the
number of subsets (typically 2,000, with in average 5 synthetic sources in
each subset), that is determinant to be used in practice. At this stage, we
still need to show that the proxy procedure we defined leads to reliable
results so that our comparisons with state-of-the-art algorithms are fair. In
particular, we aim to show that the model resulting from the proxy procedure
is not prone to over-fitting induced by an imperfect separation between
training and testing data. In the proxy procedure, over-fitting could possibly
occur if the model partially memorizes the nuisance component that is seen by
the network without the injected sources we aim to detect at inference time.
Figure 11 shows examples of detection maps obtained with the procedure of
reference described previously and its proxy version on the three datasets of
HD 95086 (2015-05-05, 2018-02-23, 2021-03-11) considered in this work. In
these experiments, we considered more than 700 synthetic sources spread over
the whole field of view. Table 6 compares quantitatively the two approaches in
terms of detection performance. These results show that the proxy procedure
leads to reliable and conservative estimations of the overall performance of
deep PACO. In addition, since synthetic sources mimicking the behaviour of
real unknown sources are recovered with comparable rates between the two
procedures, it emphasizes that the custom data-augmentation strategy of the
training step, including a random permutation of the images of the data
series, is efficient to circumvent the absence of ground-truth about real
sources. In the following, we safely use the proxy procedure to compare the
detection capability of the proposed deep PACO algorithm with other post-
processing methods.

Figure 11: Comparison between the reference procedure and its proxy version for evaluation of the detection performance of the proposed algorithm. Synthetic sources are classified as true or missed detections using a detection threshold at $\tau=0.5$. Examples of detection maps obtained in the presence of injected synthetic sources are given on three datasets of HD 95086 (2015-05-05, 2018-02-23, 2021-03-11), see Sect. 4 for the observation logs. Table 6: Comparison between the reference procedure and its proxy version for evaluation of the detection performance of the proposed algorithm. Synthetic sources are classified as true, missed, or false detections using a detection threshold at $\tau=0.5$. Mean detection results are averaged for the three datasets of HD 95086 (2015-05-05, 2018-02-23, 2021-03-11), see Sect. 4 for the observation logs. | reference procedure | proxy procedure
---|---|---
true detections | 583/728 (80.0%) | 564/728 (77.5%)
missed detections | 145/728 (20.0%) | 162/728 (22.5%)
false detections | 0 | 0

Following the previously defined proxy procedure, we present detailed results
obtained on HIP 88399 (2018-04-11), which is a SPHERE-IRDIS dataset
representative of the mean results we obtained over the eleven ones we
consider in this work. Results for the ten other datasets are reported in the
supplementary material. Figure 12 shows detection results on a sample of
10,000 synthetic sources in a diagram contrast versus angular separation for
PACO and the proposed method. Each synthetic source is classified as missed,
true, or false detection using the detection thresholds defined in Sect.
4.2.2. For PACO, setting the detection threshold at $\tau=5$ corresponds, in
average, to a realistic control (Flasseur et al., 2018a, b, c, 2020a) of the
probability of false alarms (PFA) at $5\sigma$ (i.e., PFA $\simeq 3\times
10^{-7}$). While the PFA should theoretically be controlled by the other
algorithms producing a S/N map (cADI and PCA), we have shown in previous works
(Flasseur et al., 2018a, b, c, 2020a) that the contrast curves are over-
optimistic for these algorithms (i.e., there are significantly more false
alarms than expected) due to a miss-modelling of the nuisance component. This
claim is also supported by the detection maps given in Figs. 1 to 20 of the
supplementary material. For the proposed method, converting pseudo-
probabilities into S/N scores is not feasible given that the pseudo-
probabilities are very close either to 0 or 1. For this reason, we check
empirically that the targeted false alarm rate at $5\sigma$ is satisfied. By
thresholding at $\tau=0.5$ tens of detection maps obtained with the proposed
approach, we experienced fewer false alarms than statistically expected at
$5\sigma$ so that our results can be fairly compared with the PACO results.
Figure 12 illustrates the capability of the proposed method to detect fainter
sources than PACO. From the results presented in Fig. 12, we derive the
contrast curve of each algorithm by computing, as a function of the angular
separation, the contrast value leading to an equal amount of true and missed
detections. This step is performed with a sliding window of 0.05 arcsec wide.
Figure 13 summarizes the contrast curves obtained with the five considered
algorithms. The proposed method achieves the best detection sensitivity with
an improvement in contrast up to a factor four with respect to the PACO
algorithm. We also compare the detection sensitivity of deep PACO with the
fundamental detection limit driven by the photon noise. The procedure to
compute the photon noise limit is based on a careful evaluation of the
contribution of the different sources of noise (i.e., photon, thermal
background, and detector readout noise) combined with a statistical evaluation
of the underlying S/N. This procedure will be described in details in a paper
currently in preparation. We observe that deep PACO can reach for some
datasets (see e.g. HIP 88399 (2018-04-11) in Fig. 13, and HIP 88399
(2015-05-10) as well as HIP 88399 (2016-04-16) in the supplementary material)
at large separations the best achievable detection limit driven by the photon
noise, which corresponds to an optimal unmixing between the signal of the
sources of interest and the nuisance component. Near the star, an important
gap remains (by a factor 5 to 10) between the actual performance and the
theoretical lower limit, which a sign of a lack of angular diversity in this
area. In practice, the gap between the actual contrast and the lower bound
limit depends on several characteristics such as the quality and the stability
of the observing conditions, the total amount of parallactic rotation, the
number $T$ of temporal frames, etc., as illustrated in Figs. 21 and 22 of the
supplementary materials. Reducing this gap at close separation could be
addressed by investigating ways to perform the training step of our model from
various datasets for which the presence of sources at similar locations is
very unlikely. This adaptation requires specific developments that are left
for future work.

Figure 12: Detection results for 10,000 synthetic sources in a diagram plotting contrast versus angular separation. Each synthetic source is classified as missed, true or false detection. The angular separation of false alarms is reported on the $x$-axis, and they are not associated to contrast value since they do not correspond to injected synthetic sources. Dataset: HIP 88399 (2018-04-11), see Sect. 4 for the observation logs. Figure 13: Contrast as a function of the angular separation. It can be noted that the contrast curves of cADI and PCA are over-optimistic (see text) and the angular separations leading to a PFA locally higher than ten times the targeted PFA at $5\sigma$ are not shown. The black dashed line represents the ultimate detection limit driven by the photon noise. Dataset: HIP 88399 (2018-04-11), see Sect. 4 for the observation logs. Figure 14: TPR (in percent) of deep PACO for synthetic sources missed by PACO (i.e., for $\text{S/N}\leq 5$) as a function of the angular separation (on the $x$-axis) and of the PACO’s S/N of detection (on the $y$-axis). The black dashed line represents the equivalent PACO’s detection threshold to reach $\text{TPR}=50\%$ with deep PACO (we recall that the classical detection threshold at $5\sigma$ with PACO also corresponds to $\text{TPR}=50\%$). Dataset: HIP 88399 (2018-04-11), see Sect. 4 for the observation logs. Figure 15: ROCs showing the TPR as a function of the FDR for injected synthetic sources. Dataset: HIP 88399 (2018-04-11), see Sect. 4 for the observation logs. Table 7: Mean results of AUC for ROCs giving the TPR as a function of the FDR. The scores are averaged over the eleven SPHERE-IRDIS datasets considered in this paper, see Sect. 4 for the observation logs. Only massively injected sources (10,000 per dataset) were considered. Figures 8 of (Flasseur et al., 2022) displays the corresponding ROC from which these mean results were aggregated. sep. (”) | cADI | cADI (VIP) | PCA (VIP) | PACO | proposed
---|---|---|---|---|---
$\left[0;2\right]$ | 0.53 | 0.60 | 0.62 | 0.81 | 0.91
$\left[2;4\right]$ | 0.60 | 0.65 | 0.67 | 0.83 | 0.92
$\left[4;6\right]$ | 0.11 | 0.59 | 0.55 | 0.71 | 0.88
$\left[6;7\right]$ | 0.15 | 0.38 | 0.33 | 0.77 | 0.90

Figure 14 focuses on the comparison between PACO and the proposed deep PACO
algorithm. It represents the TPR of deep PACO for synthetic sources missed by
PACO (i.e., for $\text{S/N}\leq 5$) as a function of the angular separation
and of the PACO’s S/N of detection. For instance, on this dataset about 65% of
sources below 0.5 arcsec with a S/N of detection between 3.5 and 4.0 with PACO
are detected with deep PACO (i.e., above the detection threshold $\tau=0.5$).
Similarly, more than 86% of sources between 2.0 and 4.0 arcsec with a S/N of
detection between 3.5 and 4.0 with PACO are detected with deep PACO (above the
detection threshold $\tau=0.5$). For these two examples, achieving the same
TPR with PACO would require to decrease the detection threshold at the price
to an increase of the FPR up to a mean factor of 800. Figures 23 and 24 of the
supplementary material give similar type of representation than Fig. 14 for
the ten other SPHERE-IRDIS datasets analyzed in this work.

Figure 15 gives ROCs representing the TPR as a function of the FDR for the HIP
88399 (2018-04-11) dataset. There results are obtained with the 10,000
synthetic sources considered for results presented in Fig. 12. The results are
split in four different angular separation ranges: $\left[0;2\right]$,
$\left[2;4\right]$, $\left[4;6\right]$, and $\left[6;7\right]$ arcsec. Table 7
complements this study by presenting averaged results over the eleven SPHERE-
IRDIS datasets we consider in this study, and the detailed scores for each
dataset are given in Table 2 of the supplementary material. These results
illustrate again the benefits of the proposed method in terms of precision-
recall trade-off: the AUC under ROC is improved by 9 to 17$\%$ with respect to
the best comparative algorithm for the four angular separation ranges we
consider.

####  4.2.3 Identification of candidate background point-like sources

Figure 16: Detection maps obtained with PACO and deep PACO on two epochs of
the HD 95086 star observed with SPHERE-IRDIS. Symbols classify sources as
true, missed, and false detections based on our analysis. The detection
threshold was set to $\tau=5$ for PACO producing a S/N map, and to $\tau=0.5$
for deep PACO producing a pseudo-probability map. Colors classify sources as
known real sources, candidate (background) sources, and detections with
unknown status. Datasets: HD 95086 (2015-05-05 and 2018-01-05), see Sect. 4
for the observation logs.

In this section, we take the example of a joint analysis of datasets of HD
95086 considered in this work to illustrate the ability of the proposed
approach to detect candidate faint point-like sources. Figure 16 shows the
detection maps obtained with PACO and the proposed deep PACO approach on the
2015 and 2018 observations of HD 95086. It emphasizes that two candidate
point-like sources (respectively denoted CC-11 and CC-12) are detected with
PACO at a S/N above 5 only for one of the two epochs (in 2018, which is also
one of the best SPHERE-IRDIS observations of this star) while they remain just
below the classical detection threshold at $5\sigma$ on the 2015 epoch. These
two point-like sources are detected with deep PACO on the same two epochs.
Based on their astrometry, these two faint point sources are co-moving from
one epoch to the other with the other background stars seen in the projected
field of view of the instrument, so that they could be background stars
too777The goal of this paper is not to study in details these CCs, but rather
to present a new post-processing algorithm. These CCs should be taken with
caution, and a detailed analysis is needed, in particular to exclude
systematic sources of errors.. These candidate point-like sources are not
reported in the last detailed study of the HD 95086 architecture based on the
analysis of ten SPHERE-IRDIS datasets (including the 2015 epoch but not the
2018 epoch) with classical post-processing algorithms (i.e., cADI, PCA,
TLOCI), see Fig. 2 of Chauvin et al. (2018). The analysis of more recent and
better datasets (including the 2018 one) in the SHINE survey (Langlois et al.,
2021) of the SPHERE instrument allows to identify by visual inspection (i.e.,
not based on a strict measure of S/N above the classical detection threshold
at $5\sigma$) CC-11 while CC-12 has not been identified yet. The detection of
the two CCs in the worst epoch data (2015-05-05) emphasizes again the benefits
of the proposed deep PACO algorithm. This example also illustrates the
complementarity of PACO and deep PACO. Even if deep PACO does not provide
tight confidence scores (i.e., its outputs can not be directly interpreted in
terms of a S/N), it allows to identify candidate point-like sources. Given the
locations found by deep PACO, an estimation of the PFA can be derived from the
S/N extracted on the PACO detection maps at the same locations. Using this
procedure, the theoretical PFA (not including systematic sources of errors)
for the co-localization of the candidate point sources CC-11 and CC-12 in the
two epochs is small. None of CC-11 and CC-12 are detected in the 2021 epoch of
HD 95086, as shown by Figs. 3 and 4 of the supplementary material. This
observation is consistent with the fact that the achievable contrast is worse
for this dataset than for the 2015 and 2018 observations, likely because the
observing conditions were quite average for the 2021 observations and that
only half of the total integration time was spent on the target star (the
2021’s epoch being recorded with the star-hopping technique, and the resulting
reference dataset remaining not exploited in this paper).

PACO and deep PACO detection maps from the 2015 and 2018 epochs also display a
few blobs, circled in red, above the detection threshold that are not
consistently detected on multiple epochs. These detections are located very
near the borders of the (extended) field of view. They are likely artifacts
due to the lack of temporal samples in this area to estimate the model
parameters. This hypothesis is also supported by the fact that PACO detection
maps are not stationary and does not follow a centered Gaussian distribution
with unit variance in this area. deep PACO detection map from the 2021 epoch
display a blob at 0.25”, that we likely identified as a false alarm since we
have experienced a few false alarms in the same area during training time due
to a strong stellar leakage (see Figs. 3, 4, and 22(b) of the supplementary
material).

###  4.3 Characterization results

In Sect. 4.3.1, we ground the performance of the proposed photometry
estimation module on the same datasets than the considered ones for the
detection stage by resorting to massive injections of synthetic sources. For
that purpose, we use the ARE metric (Eq. (15)) averaged per angular
separation. We also put in perspective the results of the detection and of the
characterization stages to ground the global gain brought by the proposed
algorithm. In Sect. 4.3.2, we propose a fast (and approximate) evaluation
procedure suited when large amounts of synthetic sources are considered.
Finally, in Sect. 4.3.3 we briefly discuss and compare photometry estimations
provided by our method with estimations coming from the literature on some
known sources (either exoplanets, brown dwarfs or background sources) present
in the eleven SPHERE-IRDIS datasets considered in this paper.

####  4.3.1 Characterization of synthetic sources

Figure 17: Mean ARE score on the estimated photometry of injected synthetic
sources as a function of the their contrast and of their angular separation.
From top to bottom, the panels corresponds respectively to PCA (VIP), PACO,
and the proposed algorithm. For each panel, the mean detection limit (straight
line) and the mean photon noise limit (dashed line) are superimposed. The
results are averaged azimuthally for 40,000 sources of flux drawn uniformly
between $1\times 10^{-6}$ and $3\times 10^{-5}$. Datasets: the eleven SPHERE-
IRDIS datasets considered in this work, see Sect. 4 for the recording logs.
Figure 18: Mean ARE score on the estimated photometry of injected synthetic
sources as a function of their angular separation. The results are averaged
azimuthally for sources, of flux drawn uniformly between $1\times 10^{-6}$ and
$3\times 10^{-5}$, that have been considered in Fig. 17. Note that because of
the azimuthal average, some sources belonging below the detection limit
contribute to the display results, especially at short angular separations,
see Fig. 17. Datasets: the eleven SPHERE-IRDIS datasets considered in this
work, see Sect. 4 for the recording logs.

To assess the performance of our method, we reproduce a real setting,
hereafter called reference procedure, in which we estimate the flux of
injected synthetic sources representative of (possibly unknown) exoplanets. To
do so, we first inject synthetic test sources in the considered dataset, and
then extract the patches used for training the model. It is crucial to perform
these steps in that order otherwise the model could have the opportunity to
learn the residual nuisance component below the test sources, and a bias could
be introduced.

Following this strategy, Fig. 17 shows the ARE score on the estimated
photometry of synthetic sources massively injected with the reference
procedure as a function of the their contrast and of their angular separation.
Figure 26 of the supplementary material gives the same type of plots for the
throughput (i.e., the ratio $\widehat{\alpha}/\alpha$ between the retrieved
and ground-truth source’s contrast). In both cases, results are averaged over
the eleven SPHERE-IRDIS datasets of this study. Figure 18 complements this
study by aggregating the ARE score over the source contrast. These results
show that the proposed algorithm leads, in average, to better characterization
performance than PCA (VIP) and PACO for angular separations larger than 0.8”,
with a reduction of the ARE by a factor between 1.10 and 10 with respect to
PCA (VIP), and by a factor between 1.10 and 5 with respect to PACO. Closer to
the star, the advantage is on average to PACO and to the proposed approach,
the best of the two algorithms depending on the dataset and of the source
location. The fact that PACO can perform better than the proposed algorithm is
due to the complexity to train a deep model, without leak between the train
and the test sets, from a unique dataset of interest. This effect occurs only
at short angular separations since this is the the region of the field of view
where generating multiple non-redundant training samples is the most tricky.

####  4.3.2 Efficient (and approximate) evaluation procedure

Figure 19: Comparison between the reference procedure and its proxy version
for evaluation of the characterization performance of the proposed algorithm.
(a) Mean ARE scores (see definition in Eq. (15)) are reported as a function of
the angular separation. (b) $p$-values are reported as a function of the
angular separation, and can be compared with the 5% significance threshold
given by the dashed red curve. Datasets: the eleven SPHERE-IRDIS datasets
considered in this work, see Sect. 4 for the recording logs.

The reference procedure described and applied in Sect. 4.3.1 is
computationally expensive as a new model must be trained for every test
dataset containing a dozen of synthetic test sources. As an illustration, 30
different models must be trained to get a test set of only 300 samples. In
this context and in a similar fashion to the detection stage (see Sect.
4.2.2), we also propose a fast and approximate version of the reference
procedure, that is referred as the proxy procedure in the following. Instead
of training a different model for each subset of injected sources, we train a
unique model without additional synthetic sources mimicking the behaviour of
real sources. Synthetic sources are injected a posteriori, i.e., after
training. Finally, the trained model is applied on the resulting dataset. The
gain in terms of algorithmic complexity (about a factor 2,000 with in average
5 synthetic sources in each subset) is similar to the one brought by the proxy
procedure of the detection step.

At this stage, we have to measure the ability of the proxy procedure to
provide a fair estimate of the real performance that would be achieved with
the reference procedure. Figure 19(a) compares the performance of the two
procedures in terms of mean ARE on the eleven SPHERE-IRDIS datasets considered
in this paper. In these experiments, we considered around 3,300 synthetic
sources spread over the whole field of view. We observe that the proxy
procedure leads to results very close to the ones provided by the reference
procedure, excepted near the star where the proxy procedure is over-
optimistic, i.e. a positive bias is present (up to a factor five, at worst).
This bias at short angular separations can be attributed to data leakages
between the train and the test sets in the absence of whitening procedure.
Basically, as for the detection part, the model is data-dependent and the
absence of whitening procedure as well as of the associated temporal shuffling
of the frames induces that some parts of the nuisance component are seen and
memorized by the network using the proxy procedure. This effect occurs only at
short angular separations since this is the area of the field of view where
training patches contain most likely some similar parts of the nuisance
component. This hypothesis is also supported by the absence of bias between
the proxy and the reference procedures of the detection stage. The latter
encompasses a whitening procedure and a temporal shuffling of the frames that
prevents data leakages between the train and test sets.

So, it remains to decide for a suited cut-off in terms of angular separation
to switch from the reference to the proxy procedure. For that purpose, we
resort to a binary hypothesis test. It performs a paired $t$-test (Kendall et
al., 1948), where the $\mathcal{H}_{0}$ (null) hypothesis represents the
equality between the reference and the proxy procedures while the
$\mathcal{H}_{1}$ (alternative) hypothesis represents better results with the
proxy procedure. From the results presented in Fig. 19(a), we conduct in Fig.
19(b) this statistical test on the same datasets. It displays the resulting
$p$-values as a function of the angular separation, and compares them with a
5% significance threshold. For the prescribed significance level, the two
procedures can be considered as equivalent for angular separations larger than
0.5 arcsec.

As a conclusion of this study, when the number of synthetic sources for which
the photometry should be evaluated constitutes a computational bottleneck, the
efficient and approximate proxy procedure can be safely used in the main part
of the field of view and the more expensive procedure of reference should be
used near the star.

####  4.3.3 Characterization of real sources

Table 8: Estimated photometry on some known real sources present in the considered datasets. For HD 95086, the (candidate) point-like sources are denoted as in Fig. 16. For HIP 88399 (respectively, HIP 72192), point-like sources are denoted by CC1 to CC5 (respectively, CC1 and CC2) through increasing angular separations. The astrometry (angular separation and true-North angle) has been estimated with PACO. Photometry estimations extracted from the literature has been obtained with TLOCI by (Langlois et al., 2021) for (a), TLOCI by (Desgrange et al., 2022) for (b), ANDROMEDA by (Desgrange et al., 2022) for (c), PCA by (Chauvin et al., 2017) for (e), and TLOCI by (Cheetham et al., 2019) for (f). Reported uncertainties correspond to three times the estimated standard-deviation. identification | estimated photometry $\widehat{\bm{\alpha}}$ (contrast)
---|---
target | obs. date | source | ang. sep. (”) | angle (°) | literature | PCA (VIP) | PACO | proposed
HD 95086 | 2015-05-05 | b |  0.619 $\pm$ 0.003 |  148.4 $\pm$ 0.2 |  ${}^{\text{(a)}}$(1.13 $\pm$ 0.46)$\times 10^{-5}$ |  (1.06 $\pm$ 0.06)$\times 10^{-5}$ |  (1.14 $\pm$ 0.17)$\times 10^{-5}$ |  1.11$\times 10^{-5}$
|  |  |  |  |  ${}^{\text{(b)}}$(1.10 $\pm$ 0.55)$\times 10^{-5}$ |  |  |
|  |  |  |  |  ${}^{\text{(c)}}$(1.58 $\pm$ 1.75)$\times 10^{-5}$ |  |  |
HD 95086 | 2015-05-05 | CC1 |  2.207 $\pm$ 0.003 |  35.1 $\pm$ 0.1 |  ${}^{\text{(a)}}$(6.19 $\pm$ 1.94)$\times 10^{-6}$ |  (6.33 $\pm$ 0.57)$\times 10^{-6}$ |  (7.05 $\pm$ 0.72)$\times 10^{-6}$ |  6.41$\times 10^{-6}$
HD 95086 | 2015-05-05 | CC11 |  2.346 $\pm$ 0.006 |  305.4 $\pm$ 0.2 | – |  (2.06 $\pm$ 4.61)$\times 10^{-6}$ |  (3.14 $\pm$ 0.70)$\times 10^{-6}$ |  1.80$\times 10^{-6}$
HD 95086 | 2015-05-05 | CC12 |  4.546 $\pm$ 0.010 |  52.3 $\pm$ 0.1 | – |  (2.53 $\pm$ 1.31)$\times 10^{-6}$ |  (2.70 $\pm$ 0.90)$\times 10^{-6}$ |  1.94$\times 10^{-6}$
HD 95086 | 2018-01-05 | b |  0.617 $\pm$ 0.002 |  145.7 $\pm$ 0.2 |  ${}^{\text{(b)}}$(1.20 $\pm$ 0.32)$\times 10^{-5}$ |  (1.23 $\pm$ 0.23)$\times 10^{-5}$ |  (1.20 $\pm$ 0.06)$\times 10^{-5}$ |  1.18$\times 10^{-5}$
|  |  |  |  |  ${}^{\text{(c)}}$(1.32 $\pm$ 0.35)$\times 10^{-5}$ |  |  |
HD 95086 | 2018-01-05 | CC1 |  2.241 $\pm$ 0.002 |  37.4 $\pm$ 0.1 | – |  (6.25 $\pm$ 1.88)$\times 10^{-6}$ |  (6.59 $\pm$ 0.37)$\times 10^{-6}$ |  6.39$\times 10^{-6}$
HD 95086 | 2018-01-05 | CC11 |  2.262 $\pm$ 0.004 |  305.9 $\pm$ 0.1 | – |  (2.56 $\pm$ 0.28)$\times 10^{-6}$ |  (2.75 $\pm$ 0.35)$\times 10^{-6}$ |  1.90$\times 10^{-6}$
HD 95086 | 2018-01-05 | CC12 |  4.606 $\pm$ 0.006 |  53.0 $\pm$ 0.1 | – |  (2.73 $\pm$ 3.00)$\times 10^{-6}$ |  (2.39 $\pm$ 0.40)$\times 10^{-6}$ |  2.40$\times 10^{-6}$
HD 95086 | 2021-03-11 | b |  0.630 $\pm$ 0.003 |  142.4 $\pm$ 0.2 | – |  (1.21 $\pm$ 0.23)$\times 10^{-5}$ |  (1.25 $\pm$ 0.18)$\times 10^{-5}$ |  1.20$\times 10^{-5}$
HD 95086 | 2021-03-11 | CC1 |  2.293 $\pm$ 0.005 |  40.8 $\pm$ 0.1 | – |  (5.66 $\pm$ 0.73)$\times 10^{-6}$ |  (6.30 $\pm$ 1.56)$\times 10^{-6}$ |  5.48$\times 10^{-6}$
HD 95086 | 2021-03-11 | CC11 | S/N too low | S/N too low | S/N too low | S/N too low | S/N too low | S/N too low
HD 95086 | 2021-03-11 | CC12 | S/N too low | S/N too low | S/N too low | S/N too low | S/N too low | S/N too low
HIP 88399 | 2015-05-10 | CC1 |  1.302 $\pm$ 0.002 |  210.4 $\pm$ 0.1 | – |  (7.09 $\pm$ 0.72)$\times 10^{-6}$ |  (7.39 $\pm$ 0.66)$\times 10^{-6}$ |  6.67$\times 10^{-6}$
HIP 88399 | 2015-05-10 | CC2 |  3.488 $\pm$ 0.006 |  142.4 $\pm$ 0.1 |  ${}^{\text{(a)}}$(7.94 $\pm$ 2.10)$\times 10^{-6}$ |  (7.10 $\pm$ 0.36)$\times 10^{-6}$ |  (7.56 $\pm$ 0.27)$\times 10^{-6}$ |  7.61$\times 10^{-6}$
HIP 88399 | 2015-05-10 | CC3 |  4.281 $\pm$ 0.006 |  227.7 $\pm$ 0.1 | – |  (2.36 $\pm$ 0.17)$\times 10^{-7}$ |  (3.42 $\pm$ 1.03)$\times 10^{-7}$ |  2.52$\times 10^{-7}$
HIP 88399 | 2015-05-10 | CC4 |  5.208 $\pm$ 0.003 |  322.1 $\pm$ 0.1 |  ${}^{\text{(a)}}$(5.11 $\pm$ 6.00)$\times 10^{-7}$ |  (5.57 $\pm$ 5.06)$\times 10^{-7}$ |  (4.89 $\pm$ 0.95)$\times 10^{-7}$ |  4.94$\times 10^{-7}$
HIP 88399 | 2015-05-10 | CC5 | S/N too low | S/N too low | S/N too low | S/N too low | S/N too low | S/N too low
HIP 88399 | 2016-04-16 | CC1 |  1.237 $\pm$ 0.003 |  212.8 $\pm$ 0.1 |  ${}^{\text{(a)}}$(9.38 $\pm$ 6.99)$\times 10^{-6}$ |  (1.04 $\pm$ 0.13)$\times 10^{-5}$ |  (1.22 $\pm$ 0.18)$\times 10^{-5}$ |  1.11$\times 10^{-5}$
HIP 88399 | 2016-04-16 | CC2 |  3.422 $\pm$ 0.001 |  141.7 $\pm$ 0.1 |  ${}^{\text{(a)}}$(1.32 $\pm$ 2.97)$\times 10^{-5}$ |  (1.47 $\pm$ 0.05)$\times 10^{-5}$ |  (1.64 $\pm$ 0.05)$\times 10^{-5}$ |  1.57$\times 10^{-5}$
HIP 88399 | 2016-04-16 | CC3 |  4.240 $\pm$ 0.003 |  228.6 $\pm$ 0.1 |  ${}^{\text{(a)}}$(5.35 $\pm$ 4.20)$\times 10^{-7}$ |  (5.10 $\pm$ 4.65)$\times 10^{-7}$ |  (7.64 $\pm$ 1.36)$\times 10^{-7}$ |  5.78$\times 10^{-7}$
HIP 88399 | 2016-04-16 | CC4 |  5.282 $\pm$ 0.003 |  322.6 $\pm$ 0.1 |  ${}^{\text{(a)}}$(5.11 $\pm$ 6.00)$\times 10^{-7}$ |  (7.21 $\pm$ 1.14)$\times 10^{-7}$ |  (8.02 $\pm$ 1.25)$\times 10^{-7}$ |  7.58$\times 10^{-7}$
HIP 88399 | 2016-04-16 | CC5 |  6.561 $\pm$ 0.003 |  201.8 $\pm$ 0.1 | – |  (1.10 $\pm$ 2.19)$\times 10^{-6}$ |  (1.22 $\pm$ 0.94)$\times 10^{-6}$ |  1.18$\times 10^{-6}$
HIP 88399 | 2018-04-11 | CC1 |  1.106 $\pm$ 0.001 |  217.9 $\pm$ 0.1 |  ${}^{\text{(a)}}$(1.61 $\pm$ 1.27)$\times 10^{-5}$ |  (2.01 $\pm$ 0.25)$\times 10^{-5}$ |  (2.22 $\pm$ 0.15)$\times 10^{-5}$ |  2.22$\times 10^{-5}$
HIP 88399 | 2018-04-11 | CC2 |  3.289 $\pm$ 0.001 |  140.1 $\pm$ 0.1 |  ${}^{\text{(a)}}$(2.21 $\pm$ 1.74)$\times 10^{-5}$ |  (2.65 $\pm$ 2.19)$\times 10^{-5}$ |  (2.82 $\pm$ 0.06)$\times 10^{-5}$ |  2.88$\times 10^{-5}$
HIP 88399 | 2018-04-11 | CC3 |  4.140 $\pm$ 0.003 |  230.3 $\pm$ 0.1 |  ${}^{\text{(a)}}$(8.79 $\pm$ 6.90)$\times 10^{-7}$ |  (1.10 $\pm$ 0.46)$\times 10^{-6}$ |  (9.70 $\pm$ 1.31)$\times 10^{-7}$ |  1.03$\times 10^{-6}$
HIP 88399 | 2018-04-11 | CC4 |  5.411 $\pm$ 0.002 |  323.6 $\pm$ 0.1 |  ${}^{\text{(a)}}$(1.21 $\pm$ 0.95)$\times 10^{-6}$ |  (1.31 $\pm$ 0.25)$\times 10^{-6}$ |  (1.40 $\pm$ 0.12)$\times 10^{-6}$ |  1.38$\times 10^{-6}$
HIP 88399 | 2018-04-11 | CC5 |  6.429 $\pm$ 0.002 |  202.5 $\pm$ 0.1 | – |  (4.68 $\pm$ 1.85)$\times 10^{-7}$ |  (2.83 $\pm$ 0.33)$\times 10^{-6}$ |  2.81$\times 10^{-6}$