# Deep learning for exoplanet detection and characterization by direct imaging at high contrast

Théo Bodrito, Olivier Flasseur, Julien Mairal, Jean Ponce, Maud Langlois, Anne-Marie Lagrange

► To cite this version:

Théo Bodrito, Olivier Flasseur, Julien Mairal, Jean Ponce, Maud Langlois, et al. Deep learning for exoplanet detection and characterization by direct imaging at high contrast. SF2A 2025 - Journées de la Société Française d'Astronomie & d'Astrophysique, Jul 2025, Toulouse, France. pp.1-5. hal-05281730

HAL Id: hal-05281730

https://hal.science/hal-05281730v1

Submitted on 24 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# DEEP LEARNING FOR EXOPLANET DETECTION AND CHARACTERIZATION BY DIRECT IMAGING AT HIGH CONTRAST

T. Bodrito<fnref n="1" />, O. Flasseur<fnref n="2" />, J. Mairal<fnref n="3" />, J. Ponce<fnref n="1" />,<fnref n="4" />, M. Langlois<fnref n="2" /> and A.-M. Lagrange<fnref n="5" />,<fnref n="6" />

**Abstract.** Exoplanet imaging is a major challenge in astrophysics due to the need for high angular resolution and high contrast. We present a multi-scale statistical model for the nuisance component corrupting multivariate image series at high contrast. Integrated into a learnable architecture, it leverages the physics of the problem and enables the fusion of multiple observations of the same star in a way that is optimal in terms of detection signal-to-noise ratio. Applied to data from the VLT/SPHERE instrument, the method significantly improves the detection sensitivity and the accuracy of astrometric and photometric estimation.

**Keywords:** high-contrast imaging, high angular resolution, data science, statistical methods, deep learning

## 1 Introduction

The detection of exoplanets, the characterization of their atmospheres, and the study of their formation are major challenges in modern astrophysics. Direct imaging is an observational technique of choice to address these questions but it requires to reach both a high contrast and a high angular resolution (Currie et al. 2022; Follette 2023). Alongside adaptive optics and coronagraphs, advances in data science are now critical to disentangle astrophysical signals (exoplanets and circumstellar disks) from the strong nuisance component (i.e., speckles plus noise) corrupting the observations (Pueyo 2018; Follette 2023). In this context, we present the key components of ExoMILD (Exoplanet imaging by MIxture of Learnable Distributions) introduced in our recent work (Bodrito et al. 2025a,b). ExoMILD is a deep learning approach for extracting astrophysical information from multivariate observations (spatial, temporal, spectral, multi-epoch) at high contrast. It combines statistical modeling of the nuisance with domain knowledge priors and exploits large archival database.

## 2 Proposed approach

Figure 1 summarizes the main steps of the proposed approach and can be referred to throughout this section.

### 2.1 Image formation model

An observation $y \in \mathbb{R}^{C \times T \times H \times W}$, obtained in angular and spectral differential imaging (ASDI) using the pupil tracking mode of the telescope, is composed of $T$ temporal images of $H \times W$ pixels in $C$ spectral channels. The forward model is:

$$ y_{c,t} = \beta_c \mathbf{D}_{\lambda_c/\lambda_1} s_{1,t} + \sum_{k=1}^{K} \alpha_c^{(k)} h_c(x_t^{(k)}) + \epsilon_{c,t}, \quad (2.1) $$

where $\beta_c \in \mathbb{R}$ is the speckle amplitude in channel $c$, $\mathbf{D}_{\lambda_c/\lambda_1}$ is the operator that aligns speckles from the reference channel ($\lambda_1$) to those of channel $c$ via a spatial homothety with factor $\lambda_c/\lambda_1$, $s_{1,t} \in \mathbb{R}^{H \times W}$ is the

<footnote n="1">¹ Département d’Informatique de l’École normale supérieure (ENS-PSL, CNRS, Inria)</footnote>

<footnote n="2">² Université Claude Bernard Lyon 1, Centre de Recherche Astrophysique de Lyon UMR 5574, ENS de Lyon, CNRS, Villeurbanne, France</footnote>

<footnote n="3">³ Université Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK</footnote>

<footnote n="4">⁴ Courant Institute and Center for Data Science, New York University</footnote>

<footnote n="5">⁵ LIRA, Observatoire de Paris, Université PSL, Sorbonne Université, Université Paris Cité, CNRS, Meudon, France</footnote>

<footnote n="6">⁶ Université Grenoble Alpes, Institut de Planétologie et d'Astrophysique de Grenoble</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Proposed approach combining estimated and learnable parameters. It exploits the spectral diversity of the observations $\mathbf{y}$ through a convolutional mixture model that captures nuisance symmetries and correlations at different spatial scales. The network $f_\nu$ (U-Net) learns a prior on the exoplanet signals to filter the detection score $\hat{\gamma}$ produced by the statistical model.</figcaption></figure>

speckle contribution at $\lambda_1$, and $\epsilon_{c,t} \in \mathbb{R}^{H \times W}$ is an additive noise term (thermal, photon, readout), with $\epsilon \ll s$ near the star. The number of exoplanets $K$ is unknown. The spatial signature of exoplanet $k$ in channel $c$ is given by the instrumental PSF $h_c \in \mathbb{R}^{H \times W}$ with flux $\alpha_c^{(k)} \in \mathbb{R}_+$. Its position $\mathbf{x}_t^{(k)}$ at time $t$ is defined by an apparent rotation $r(\mathbf{x}_1^{(k)}, \phi_t)$ around the star (induced by ASDI), where $\phi_t$ is the cumulative parallactic rotation angle between $t_1$ and $t$, and $\mathbf{x}_1^{(k)} \in \mathbb{R}^2$ is its initial position in the image at $t=1$.

## 2.2 Convolutional statistical model

We propose a local representation of the nuisance component. Rather than modeling pixel correlations directly at the patch level, we represent them using a multivariate Gaussian model in a feature space obtained through learnable linear projection:

$$ \forall (j, c, t) \in [1, M] \times [1, L] \times [1, T], \quad \mathbf{A} \mathbf{y}_{c,t}^{(j)} \sim \mathcal{N}(\mathbf{m}_{c,j}, \mathbf{C}_j), $$

where $\mathbf{y}_{c,t}^{(j)} \in \mathbb{R}^p$ is a patch $j$ (among the $M$ possible) extracted from the observations $\mathbf{y}_{c,t}$, and $\mathbf{A} \in \mathbb{R}^{m \times p}$ is a learnable projection matrix with $m \le p$. The parameters $\mathbf{m}_{c,j} \in \mathbb{R}^p$ and $\mathbf{C}_j \in \mathbb{R}^{p \times p}$, representing respectively the temporal mean in channel $c$ and the spatial covariance of the Gaussian distribution, are estimated by maximum likelihood with an unsupervised shrinkage regularization of $\mathbf{C}_j$, ensuring an optimal bias/variance trade-off (Flasseur et al. 2024). We denote $\mathbf{\Omega}_{c,j} = \{\mathbf{m}_{c,j}, \mathbf{C}_j\}$. The linear projection $\mathbf{A}$, learned from the data, decorrelates the feature space where the statistical distribution is defined from the pixel space. This model captures spatial correlations at longer spatial ranges than those allowed by the maximum patch size in the pixel space, thereby maintaining reliable covariance estimation. Assuming independence of patch collections along an exoplanet trajectory, the local statistical model provides an estimate of the time-invariant flux $\hat{\alpha}_c$ at the initial position $\mathbf{x}_1 \in \mathbb{R}^2$ by maximizing the global likelihood:

$$ \hat{\alpha}_c = \arg \max_{\alpha_c} \ell(\alpha_c, \mathbf{x}_1), \quad \text{with} \quad \ell(\alpha_c, \mathbf{x}_1) = \prod_{c,t} \prod_{j \in S(\mathbf{x}_t)} \mathbb{P} \left( \mathbf{A} \left( \mathbf{y}_{c,t}^{(j)} - \alpha_c \mathbf{h}_c^{(j)}(\mathbf{x}_t) \right) \middle| \hat{\mathbf{\Omega}}_{c,j} \right)^{w_j} \quad (2.2) $$

where $\mathbf{h}^{(j)}(\mathbf{x}_t) \in \mathbb{R}^p$ is the patch $j$ of the PSF centered at $\mathbf{x}_t$, and $S(\mathbf{x}_t)$ is the subset of modeled distributions for patches around position $\mathbf{x}_t$, each distribution $j$ being weighted by $w_j \in \mathbb{R}_+$ with $\sum_j w_j = 1$. This approach is robust because it aggregates contributions from all patches partially overlapping at location $\mathbf{x}_t$, thus modeling the nuisance as a convolutional process influenced by multiple overlapping noise sources. Solving (2.2) at position $\mathbf{x}_1$ yields the flux estimator $\hat{\alpha}_c$, its associated standard deviation $\hat{\sigma}_c$, and the detection criterion $\hat{\gamma}$, which allows to decide between the absence (hypothesis $\mathcal{H}_0$) or presence (hypothesis $\mathcal{H}_1$) of an exoplanet located at $\mathbf{x}_1$ in the image at $t=1$ by comparison with a detection threshold $\tau$:

$$ \hat{\alpha}_c = \frac{\sum_{c,t} b_{c,t}(\mathbf{x}_t)}{\sum_{c,t} a_{c,t}(\mathbf{x}_t)}, \quad \hat{\sigma}_c = \frac{1}{\sqrt{\sum_{c,t} a_{c,t}(\mathbf{x}_t)}}, \quad \hat{\gamma} = \frac{\sum_{c,t} b_{c,t}(\mathbf{x}_t)}{\sqrt{\sum_{c,t} a_{c,t}(\mathbf{x}_t)}}, \quad \mathcal{H}_1 \stackrel{\tau}{\ge} \mathcal{H}_0 \quad (2.3) $$
with

$$b_{c,t}(\mathbf{x}_t) = \sum_{j \in S(\mathbf{x}_t)} w_j \mathbf{h}_c^{(j)}(\mathbf{x}_t)^\top \mathbf{A}_j^\top \hat{\mathbf{C}}_j^{-1} (\mathbf{A}_j \mathbf{y}_{c,t}^{(j)} - \hat{\mathbf{m}}_{c,j}), \quad a_{c,t}(\mathbf{x}_t) = \sum_{j \in S(\mathbf{x}_t)} w_j \mathbf{h}_c^{(j)}(\mathbf{x}_t)^\top \mathbf{A}_j^\top \hat{\mathbf{C}}_j^{-1} \mathbf{A}_j \mathbf{h}_c^{(j)}(\mathbf{x}_t).$$

The detection criterion $\hat{\gamma}$ corresponds to the generalized likelihood ratio test (GLRT) and follows, under $H_0$, a Gaussian distribution $N(0, 1)$. It can be interpreted as a detection signal-to-noise ratio. In practice, the estimates of $\hat{\alpha}_c$, $\hat{\sigma}_c$, and $\hat{\gamma}$ are jointly performed with those of the statistics $\hat{\Omega}_{c,j}$ characterizing the nuisance, since the latter may initially be biased by the presence of an exoplanet. The estimators $\hat{\alpha}_c$, $\hat{\sigma}_c$, and $\hat{\gamma}$ in $\mathbb{R}^{H \times W}$ are obtained at every point in the field by solving (2.2) independently for each location of the field of view. After thresholding the detection map $\hat{\gamma}$ at a false alarm level $1 - \Phi(\tau)$, with $\Phi$ denoting the Normal cumulative distribution function, a refined estimation of the flux and subpixel position of exoplanets is obtained by iterative optimization around the maximum of the GLRT.

## 2.3 Mixture of distributions

**Multiple scales modeling.** We enrich the model introduced in Sect. 2.2 by integrating distributions at multiple spatial scales using patches of variable sizes $p \in \mathcal{P}$. The set of modeled distributions is denoted $\cup_{p \in \mathcal{P}} S_p(\mathbf{x}_t)$, where $S_p(\mathbf{x}_t)$ gathers the patches of size $p$ that contain $\mathbf{x}_t$. The linear projection $\mathbf{A}$ reduces the feature dimension for large patches ($m < p$), ensuring the scalability of the method.

**Symmetries modeling.** In direct imaging, speckles display central and rotational symmetries (see Fig. 1, left) that arise from instrument physics. We leverage these symmetries to model the nuisance and to mitigate the effect of self-subtraction caused by the limited apparent rotation of exoplanets near the star. This limited diversity leads to a bias in estimated parameters $\hat{\Omega}_{c,j}$ due to the possible presence of exoplanets at position $j$, which in turn reduces detection sensitivity. To overcome this issue, we model the joint distribution of patches harboring similar speckle structures extracted at $j$ after rotating the images by $2\pi/K$:

$$\forall (j, c, t), \quad \mathbf{A}[\mathbf{R}_{2\pi k/K}(\mathbf{y}_{c,t})^{(j)}]_{k=1:K} \sim \mathcal{N}(\mathbf{m}_{c,j}, \mathbf{C}_j), \qquad (2.4)$$

with $\mathbf{A} \in \mathbb{R}^{m \times Kp}$. The parameters of the resulting mixture model (for $K = 1, 2, 4$) are less sensitive to self-subtraction, since the probability of observing an exoplanet simultaneously in all modeled patches is low.

**Joint spatio-spectral modeling.** Since the nuisance in any channel $c$ is related to that in channel $c = 1$ (Sect. 2.1), we propose a joint spatio-spectral model:

$$\forall (j, c, t), \quad \mathbf{A} \beta_{c,j}^{-1} \mathbf{D}_{\lambda_1/\lambda_c}(\mathbf{y}_{c,t})^{(j)} \sim \mathcal{N}(\mathbf{m}_{c,j}, \mathbf{C}_j). \qquad (2.5)$$

This approach reduces the bias on the estimation of the parameters $\hat{\Omega}_{c,j}$ by distinguishing speckles from exoplanets through their different spectral signatures. The local amplitude $\hat{\beta}_{c,j}$ is estimated via the standard deviation of the pixels in the collection $\{\mathbf{D}_{\lambda_1/\lambda_c}(\mathbf{y}_{c,t})^{(j)}\}_t$.

## 2.4 Supervised learning

We combine the statistical nuisance model (Sects. 2.3–2.4) with a prior knowledge on exoplanet signals. The spatial distribution of flux $\hat{\alpha}_c$ is thus estimated as:

$$\hat{\alpha}_c = \arg \min_{\alpha_c \in \mathbb{R}^{H \times W}} \underbrace{\varphi_{\theta}(\alpha_c, \chi(\mathbf{y}_c))}_{\text{data fidelity}} + \underbrace{\psi_{\nu}(\alpha_c)}_{\text{exoplanet prior}}. \qquad (2.6)$$

The data fidelity term $\varphi_{\theta}$ is defined by the joint log-likelihood of the nuisance model (Sects. 2.2–2.3), and $\psi_{\nu}$ is a prior on exoplanet signals implemented by a neural network. The parameters $\theta = \{\mathbf{A}_j, w_j\}_j$ and $\nu$ are jointly learnt from training data, while the nuisance parameters $\chi(\mathbf{y}_c) = \{\hat{\Omega}_{c,j}, \hat{\beta}_{c,j}\}$ are specific to each observation $\mathbf{y}_c$. We propose a two-step approach to solve (2.6) and estimate $\{\hat{\alpha}_c, \hat{\sigma}_c, \hat{\gamma}\}$. The first step solves it using only the data fidelity term, and the second step applies the neural network alone, thereby filtering the flux estimated from the previous step. This two-step approach relies on the fact that, under $H_0$, the criterion $\hat{\gamma}$ approximately follows a Gaussian distribution $N(0, 1)$. Thus, extracting exoplanet signals amounts to denoising $\hat{\gamma}$ in order to remove residual noise. This is formalized as $\tilde{\gamma} = f_{\nu}(\hat{\gamma})$, where $f_{\nu}$ is a denoiser implemented by a U-Net and $\tilde{\gamma}$
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Examples of detection maps $\tilde{\gamma}$ obtained by injecting synthetic exoplanets into real data. Arrows indicate detections achieved only with the proposed approach.</figcaption></figure>

<table><caption>**Table 1.** Aggregated detection scores (AUC) for different algorithms in ADI (independent single-spectral processing) and ASDI (joint multi-spectral processing). The ASDI mode is not supported by MODEL&CO.</caption><thead><tr><th>algorithm</th><th>mono-spectral (ADI)</th><th>multi-spectral (ASDI)</th></tr></thead><tbody><tr><td>cADI</td><td>0.369 ± 0.010</td><td>0.416 ± 0.011</td></tr><tr><td>KLIP</td><td>0.478 ± 0.010</td><td>0.515 ± 0.009</td></tr><tr><td>PACO</td><td>0.520 ± 0.011</td><td>0.693 ± 0.011</td></tr><tr><td>MODEL&amp;CO</td><td>0.644 ± 0.010</td><td>N/A</td></tr><tr><td><strong>proposed</strong></td><td><strong>0.645 ± 0.005</strong></td><td><strong>0.762 ± 0.008</strong></td></tr></tbody></table>

is the final detection map. The flux uncertainty is assumed to be only marginally affected by the network (i.e., $\tilde{\sigma}_c = \hat{\sigma}_c$), since it is mainly driven by speckle variability, which is already well captured by the statistical model in the first step. The estimated spatial distribution of flux is then given by $\tilde{\alpha}_c = f_{\nu}(\tilde{\gamma}) \times \hat{\sigma}_c$. This formulation yields a pixel-wise Gaussian distribution: $N(\tilde{\alpha}_c, \tilde{\sigma}_c)$.

The parameters $\theta$ and $\nu$ are estimated by supervised learning using realistically simulated synthetic exoplanets in real data via (2.1), with a ground-truth flux distribution $\alpha^{gt}$. The objective function minimized is the joint log-likelihood between $\alpha^{gt}$ and the estimates $(\tilde{\alpha}, \tilde{\sigma})$:

$$ \mathcal{L}(\tilde{\alpha}, \tilde{\sigma}, \alpha_{gt}) = 0.5 \sum_c (\tilde{\alpha}_c - \alpha_c^{gt})^2 / \tilde{\sigma}_c^2 + \log \tilde{\sigma}_c. \quad (2.7) $$

Finally, the model is calibrated to associate $\tilde{\gamma}$ with a false alarm probability. This is done by estimating the cumulative distribution function of $\tilde{\gamma}$ under the null hypothesis $H_0$ from a non-overlapping calibration dataset.

## 2.5 Optimal fusion of multiple observations

To further improve detection sensitivity, we propose to fuse $N$ observations $\mathbf{y}$ of the same star acquired over several years and processed with ExoMILD. While the orbital motion of exoplanets is negligible over the few hours required for a single observation $\mathbf{y}$, it becomes significant over several months. Multi-observation fusion therefore amounts to jointly estimating the orbital parameters $\mu \in \mathbb{R}^7$ of an exoplanet together with its detection (Dallant et al. 2023). The optimal fusion of the $N$ observations $\{\mathbf{y}_n\}_{n=1:N}$ amounts to combining the sufficient statistics $\{a_n, b_n\}_{n=1:N}$ produced for each observation. The maximum likelihood estimator of the orbital elements $\mu$ is:

$$ \hat{\mu} = \arg\max_{\mu} \left\{ C(\mu) = \sum_{n=1}^{N} \frac{\left( [\mathbf{b}_n(\mathbf{x}_t(\mu))]_+ \right)^2}{a_n(\mathbf{x}_t(\mu))} \right\}. \quad (2.8) $$

This estimator enables the optimal combination of the signal of a potential exoplanet along its orbit across the $N$ observations. The fused detection criterion, corresponding to the multi-observation GLRT, can be approximated as $\gamma^{comb}(\mu) \approx \sqrt{C(\mu)}$. The resolution of (2.8) by exploring the orbital parameters $\mu$ can be performed using a Hamiltonian Monte Carlo algorithm, as in the PACOME algorithm (Dallant et al. 2023).

## 3 Results

We evaluate the performance of the method on data from the VLT/SPHERE instrument. The 4D data $\mathbf{y}$ consist of $T \in [15; 300]$ images of $H \times W = 256^2$ pixels in $L = 2$ spectral channels. The model is trained on 220 observations from the instrument's archives. Detection performance is assessed using the area under the receiver operating characteristic curve (AUC), which quantifies the precision-recall trade-off as a function of the detection threshold. We compare the proposed approach with cADI (Marois et al. 2006), KLIP (Soummer
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Fusion of 23 observations of the star HR 8799. Left: Examples of individual detection maps $\hat{\gamma}$ and fused maps $\hat{\gamma}_{comb}$ centered on the three known exoplanets. Right: Orbits $\hat{\mu}$ estimated jointly with the detection.</figcaption></figure>

et al. 2012; Amara & Quanz 2012), PACO (Flasseur et al. 2020a,b), and MODEL&CO (Bodrito et al. 2024). cADI and KLIP empirically estimate the nuisance and then subtract it from the data. cADI uses a temporal median (since speckles are temporally quasi-static), whereas KLIP applies a truncated principal component decomposition. PACO statistically models the spatial correlations of the nuisance, capturing its local structure but ignoring large-scale correlations. All these approaches remain limited by self-subtraction near the star, which reduces detection sensitivity. MODEL&CO addresses this limitation with a nuisance model learned from multiple observations via deep learning. Our approach follows this line but provides better interpretability through statistical modeling, explicitly incorporating symmetries, multi-scale correlations, and spectral diversity.

Figure 2 displays detection maps on SPHERE data with synthetic exoplanets simulated using (2.1) for increased diversity. Table 1 reports the averaged detection metrics over 5 datasets. The proposed ExoMILD method significantly improves performance, with an additional gain provided by joint multi-spectral processing. Finally, Fig. 3 illustrates the fusion strategy described in Sect. 2.5 by combining 23 observations of the star HR 8799 collected between 2014 and 2022. The proposed approach enables the detection of the three known exoplanets in the field of view and the joint estimation of their orbits, while also improving detection confidence.

## 4 Conclusions

The proposed approach is tailored to the specific challenges of direct imaging: (i) very low signal-to-noise ratios and non-stationary noise, (ii) detection of rare events, and (iii) absence of ground truth. Using real data from the VLT/SPHERE instrument, we show that this approach enables fine modeling and effective subtraction of the nuisance component, leading to reliable and nearly optimal estimates of the astrophysical quantities of interest. This results in significantly improved detection sensitivity and more accurate astro-photometric characterization. The proposed approach is also scalable and readily applicable to large-scale surveys.

This work was supported by the French government under the France 2030 program (PR[AI]RIE-PSAI, ANR-23-IACL-0008; MIAI 3IA Institute, ANR-19-P3IA-0003; PEPR Origins, ANR-22-EXOR-0016), the European Research Council (ERC grants APHELEIA 101087696 and COBREX 885593), the French National Programs (PNP, PNPS), and the CNRS/INSU Action Spécifique Haute Résolution Angulaire (ASHRA) co-funded by CNES.

## References

* Amara, A. & Quanz, S. P. 2012, Monthly Notices of the Royal Astronomical Society, 427, 948

* Bodrito, T., Flasseur, O., Mairal, J., et al. 2024, Monthly Notices of the Royal Astronomical Society, 534, 1569

* Bodrito, T., Flasseur, O., Mairal, J., et al. 2025a, in XXXe Colloque GRETSI

* Bodrito, T., Flasseur, O., Mairal, J., et al. 2025b, in Computer Vision and Pattern Recognition Conference

* Currie, T., Biller, B., Lagrange, A.-M., et al. 2022, arXiv preprint arXiv:2205.05696

* Dallant, J., Langlois, M., Flasseur, O., et al. 2023, Astronomy & Astrophysics, 679, A38

* Flasseur, O., Denis, L., Thiébaut, É., et al. 2020a, Astronomy & Astrophysics, 637, A9

* Flasseur, O., Denis, L., Thiébaut, É., et al. 2020b, Astronomy & Astrophysics, 634, A2

* Flasseur, O., Thiébaut, E., Denis, L., et al. 2024, in 32nd European Signal Processing Conference (EUSIPCO), IEEE

* Follette, K. B. 2023, Publications of the Astronomical Society of the Pacific, 135, 093001

* Marois, C., Lafrenière, D., Doyon, R., et al. 2006, The Astrophysical Journal, 641, 556

* Pueyo, L. 2018, Handbook of Exoplanets, 705

* Soummer, R., Pueyo, L., & Larkin, J. 2012, The Astrophysical Journal Letters, 755, L28