Direct imaging: subtracting a star to see its planets

In October 2007 Christian Marois and his collaborators were looking at Gemini North images of the young star HR 8799. Two faint points sat near the star. In any single exposure they were no more convincing than dozens of speckles, the blobs of starlight that a telescope’s residual optical errors scatter around every bright star. What made them different was how they moved. Over the night the sky turned relative to the instrument, so the two points rotated around the star while the speckles, fixed to the optics, stayed put. A third point appeared in Keck data in 2008. Going back through the archive, the team then found planets b and c in Keck images from 2004, taken without the observing technique that made them visible in 2007. The 2008 paper reported the first directly imaged multi-planet system (Marois et al. 2008).

The story contains most of this chapter. The planets were there in old data that had not been processed in the right way. They were found by their motion relative to the noise, not by their brightness. Their reality was settled by later epochs that showed common proper motion with the star and orbital motion around it. Seeing, finding and confirming happened at different times.

Direct imaging means measuring photons from the planet separately from its star. That sounds simpler than inferring a planet from a stellar wobble or a transit, and in one sense it is: the planet’s position, colour and spectrum are measured directly. The difficulty is contrast. A young giant planet glowing in the near infrared may be $10^{-5}$ to $10^{-6}$ of its star’s brightness. A mature giant seen in reflected light is fainter still. The star’s light is spread by diffraction and by imperfect optics over the very region where planets live. This chapter is about how that light is removed, how the remainder is tested for planets, and how the test is extended across several nights by following the planet along its orbit.

The problem: a planet under the speckles

Diffraction, speckles and contrast

A perfect telescope of diameter $D$ images a star as an Airy pattern with a core of width about $\lambda/D$ and rings that fall off slowly. Adaptive optics corrects most of the atmospheric wavefront error, but not all. What remains, together with small static errors in the telescope and instrument, scatters starlight into a field of speckles. Each speckle has roughly the size and shape of the instrument’s point-spread function (PSF), so it looks exactly like a faint companion.

Optical suppression comes first. Bernard Lyot invented the coronagraph to observe the solar corona without an eclipse, and installed his first instrument at Pic du Midi on 25 July 1930. His 1939 George Darwin Lecture explains why simply blocking the Sun’s disk failed: diffraction at the edge of the objective and scattering in the optics still flooded the field. His answer was an occulting element followed by a pupil stop that removes the diffracted light (Lyot 1939). Stellar coronagraphs keep this division of labour, with apodizers, phase masks and vortex masks in place of a simple disk. None of them removes changing aberrations, centring errors or the photon noise of leaked light. Numerical PSF subtraction became essential because the residual starlight was neither stable nor known in advance.

“Contrast” is used for at least four different quantities, and confusing them is the commonest error in comparisons between methods. It can mean the optical attenuation of the coronagraph; the standard deviation of residuals after subtraction, divided by the stellar flux; a throughput-corrected detection threshold at some nominal significance; or the flux at which injected planets are recovered with a stated completeness and a stated false-alarm rate. The original angular differential imaging paper reached a difference of about 20 magnitudes from Vega beyond 8 arcseconds (Marois et al. 2006). The final analysis of the SPHERE SHINE survey reports median contrasts near $10^{-6}$ at 0.8 arcsec and $2\times10^{-7}$ beyond 2.5 arcsec (Chomez et al. 2025). These are different instruments, wavelengths and threshold conventions, not points on one technology curve.

The signal model

Let $\mathcal C_\lambda$ denote propagation through telescope and coronagraph, $P(\mathbf u)$ the pupil, and $a_t$, $\phi_t$ amplitude and phase errors at time $t$. The stellar field in the focal plane is

and a useful model of the mean detected intensity is

Here $F_\star$ is the stellar flux, $f_\lambda$ the planet’s flux, $h$ its off-axis response (including coronagraph transmission and detector sampling) at sky position $\boldsymbol\rho$, and $B$ the background. Star and planet are mutually incoherent, so their intensities add. For small aberrations, writing $E_\star=E_0+\delta E$,

A small change in the wavefront therefore changes the intensity coherently over many pixels. This is the physical reason a low-dimensional basis can describe the speckle pattern, and the reason principal-component methods work at all. It is a motivation, not a guarantee of low rank.

Two facts about noise govern everything that follows. First, subtracting a perfect model of the mean stellar pattern does not remove its photon noise: the variance of the leaked starlight stays in the residual. Second, differencing two noisy exposures adds their variances. A reference image reduces our uncertainty about the stellar pattern but contributes its own noise. No algorithm restores photons the optics threw away.

Diversity: the planet must move differently from the speckles

Subtraction needs some property that distinguishes planet from star. The observation supplies it; the algorithm only exploits it. Three kinds of diversity are used.

Angular differential imaging (ADI) keeps the telescope pupil fixed on the detector and lets the sky rotate. Speckles tied to the optics stay put, while a planet at angular radius $\rho$ moves along an arc. After a rotation $\Delta\theta$ it has moved

A reference frame is safe to use only if the planet has moved by more than some protection distance $qw$, where $w$ is the PSF width. With $q=1$ and $w=\lambda/D$, a planet at $10\lambda/D$ needs about $0.1$ radian, or $5.7^\circ$, of rotation, but a planet at $2\lambda/D$ needs half a radian, nearly $29^\circ$ (an illustrative calculation). The required rotation rises steeply inward, which is why the innermost separations are the hardest. The difficulty is partly informational: with little rotation, a planet and a slowly drifting speckle become nearly the same signal.

Spectral differential imaging (SDI) uses wavelength. Speckles produced by the pupil scale with wavelength, $\rho_{\rm sp}\propto\lambda$, while a planet stays at a fixed sky position. Magnifying each wavelength image by $\lambda_0/\lambda$ aligns the speckles and moves the planet to $\boldsymbol\rho\lambda_0/\lambda$, a displacement

For wavelength ratios $\lambda/\lambda_0=1$, 1.1 and 1.2, a planet at radius 5 (in any unit) appears at 5, 4.55 and 4.17 after speckle alignment. Methane-band imaging was attractive because cool planets absorb strongly where their stars do not, so one channel can serve as a nearly planet-free reference (Flasseur et al. 2020).

Reference differential imaging (RDI) uses other stars, observed under similar conditions, as the source of the stellar pattern. The target’s planet cannot contaminate the reference, but no other star matches perfectly in colour, centring, wavefront and thermal state.

One statistic underneath many algorithms

Raw speckle intensity is not Gaussian. A common model is the modified Rician distribution,

with $I_c$ the coherent and $I_s$ the speckle intensity, mean $I_c+I_s$ and variance $I_s^2+2I_cI_s$; $I_0$ is a modified Bessel function. Its long tail is why bright speckles look like planets. Subtraction and averaging can make residuals more nearly Gaussian, but that has to be checked for each reduction (Mawet et al. 2014).

Once residuals are modelled as Gaussian with covariance $C$, the detection of a planet with known processed template $\mathbf p$ in residual $\mathbf y=f\mathbf p+\mathbf n$ reduces to two numbers:

For a non-negative flux, twice the log-likelihood ratio is $\max(0,S)^2$. Independent exposures add their signed $A$ and their positive $B$. This is the same matched filter that Chapter 1 derives and that every stacking method in Part V uses. Matched filtering of a residual image, ANDROMEDA, PACO and the orbit-stacking engines below all compute $A$ and $B$; they differ in how they estimate $C$ and $\mathbf p$.

The template matters as much as the covariance. If a linear processing operator $L$ is applied to the data, the planet and the noise are transformed together: $\mathbf p_{\rm out}=L\mathbf p_{\rm in}$ and $C_{\rm out}=LC_{\rm in}L^T$. Filtering a PSF-subtracted image with an unprocessed, positive PSF ignores the negative lobes that subtraction carves around a planet. Treating resampled pixels as independent overcounts information.

Small samples: why five sigma is not five sigma close to the star

At a radius of $r$ resolution elements ($r$ in units of $\lambda/D$), an annulus contains only about $n=2\pi r$ independent speckle realisations. More detector pixels do not add more independent speckles. Dimitri Mawet and collaborators pointed out in 2014 that the usual procedure, comparing a candidate aperture with the standard deviation of other apertures at the same radius, is then a small-sample test. Comparing a candidate value $x_1$ with the mean $\bar x_2$ and sample standard deviation $s_2$ of $n_2=n-1$ background elements gives

which, for independent Gaussian samples with a common unknown variance, follows Student’s distribution with $n_2-1$ degrees of freedom (Mawet et al. 2014). The heavier tail and the extra factor both raise the threshold.

The effect is not subtle. At the false-positive fraction of a one-sided Gaussian $5\sigma$ test, $2.87\times10^{-7}$ per test, an annulus with $n=6$ elements (about $1\lambda/D$) requires a candidate 62 times the measured background standard deviation; $n=12$ (about $2\lambda/D$) requires 11.7; only by $n\approx100$ does the threshold fall to 5.4. Many published "$5\sigma$" contrast curves at small separations were therefore far less conservative than their label. The correction handles one specific problem, independent Gaussian elements with unknown variance. It does not repair heavy tails, noise models fitted to the same data, searches over many spectral templates, or a person choosing the reduction that makes a candidate look best.

The other ingredients of a contrast curve are throughput and completeness. A throughput-corrected threshold has the form

with $\tau$ the threshold, $\sigma_{\rm aperture}$ the noise in a planet-sized aperture, and $T_{\rm alg}$, $T_{\rm coro}$ the algorithmic and coronagraphic throughputs. A planet whose mean signal sits exactly at the threshold is recovered half the time: the threshold is a 50% completeness level, not a sensitivity floor. And the false-alarm probability of a survey is not that of a single test. For $M$ independent tests, $P_{\rm FA}=1-(1-\alpha)^M\simeq M\alpha$, and positions, spectral templates, numbers of retained modes and inspected reductions all count.

Subtracting the star

Coronagraphs and calibration are part of the model

A planet near the edge of a coronagraphic mask or a phase transition has a position-dependent and often asymmetric response. The 2025 JWST/MIRI observation of TWA 7 shows how much this matters. Its two telescope rolls placed the candidate at different distances from a transition of the four-quadrant phase mask, giving modelled transmissions of about 0.62–0.66 in one roll and 0.28–0.31 in the other; the signal-to-noise ratios in the two reductions were 30 and 18, and the authors used the favourable roll for photometry (Lagrange et al. 2025). Equal exposure time did not mean equal information. A reproducible reduction records the unocculted PSF or satellite-spot calibration, centring, wavelength solution, parallactic angles and masks, and decides whether to resample the data or carry the model onto the original detector grid.

Classical ADI

Marois, David Lafrenière, René Doyon, Bruce Macintosh and Daniel Nadeau published angular differential imaging in 2006, from Gemini North’s Altair adaptive optics and the NIRI camera. The target became its own PSF reference, avoiding the colour and brightness mismatch of a separate star. Single subtractions attenuated quasi-static structure by about a factor of five, and rotated coaddition suppressed it further (Marois et al. 2006). Their recipe was specific: remove an all-frame median, work in 0.6-arcsec annuli, choose four suitably separated frames per annulus, median-combine them into a reference, subtract, derotate and median-combine.

calibrate and register the cube
for each target frame t:
    for each annulus a:
        J = frames whose parallactic angle moves a planet at radius(a) by > q * FWHM
        reference = median of the frames in J nearest in time to t
        residual[t, a] = data[t, a] - fitted_scale * reference
derotate residuals by their parallactic angles
combine in sky coordinates; calibrate planet throughput by injection

Median construction costs $O(TP)$ for $T$ frames of $P$ pixels with linear-time selection, and derotation another $O(TP)$. The design tension is visible in the selection rule: frames far apart in time give good protection to the planet but poorly correlated speckles; frames close in time match the speckles but subtract part of the planet. ADI also attacks extended disks, which occupy the same detector pixels in every frame, producing narrow artificial features and negative shadows. A clean-looking residual image is not proof of a faithful one.

LOCI: a local least-squares reference

In 2007 Lafrenière and collaborators replaced the median by the locally optimal linear combination of reference images (LOCI). Let $X$ hold $K$ reference images over an optimisation zone of $P_o$ pixels, $\mathbf y$ the target, and $W$ a weight or mask matrix. Solve

and subtract $X_s\widehat{\mathbf c}$ only within a smaller subtraction zone. The large fitting zone constrains the speckles while diluting the influence of a planet inside it. Their demonstration used optimisation zones of $N_A=300$ PSF cores and a minimum displacement parameter $N_\delta=0.5$; it reached throughput-corrected limits of 13.9, 16.1 and 16.9 magnitudes at 1, 2 and 3 arcsec, with speckle attenuation of about 10–12 per subtraction, and up to three times better sensitivity at small separation than the preceding ADI algorithm (Lafrenière et al. 2007). Those limits used the historical $5\sigma$ convention; Mawet’s small-sample argument came later.

Per zone, forming the Gram matrix costs $O(P_oK^2)$, a dense solve $O(K^3)$. The normal equations square the condition number, so QR, SVD or a regularised solve is preferable to an explicit inverse. ADI with optimised subtraction of this kind is what produced the HR 8799 discovery.

KLIP and PCA

Two independent 2012 papers recast the reference problem as principal-component analysis. Rémi Soummer, Laurent Pueyo and James Larkin’s Karhunen–Loève Image Projection (KLIP) was demonstrated on archival HST/NICMOS data (Soummer et al. 2012). Adam Amara and Sascha Quanz’s PynPoint was demonstrated on VLT/NACO data, with a gain of up to a factor of five near 0.29 arcsec relative to LOCI in its injection test (Amara & Quanz 2012).

For mean-subtracted references $R\in\mathbb R^{K\times P_z}$ in a zone of $P_z$ pixels, solve the small eigenproblem and form normalised eigenimages,

keep the first $q$ of them as the columns of $Z_q$, and project them out of the target:

The ordering by eigenvalue minimises the expected squared reconstruction error of the reference set. It does not maximise planet S/N, and minimum residual variance is not the goal. For a fixed basis the planet’s response is exact, $(I-Z_qZ_q^T)\mathbf p$, which is the starting point of forward modelling.

The cost is dominated by the covariance, $O(P_zK^2)$, the eigendecomposition, $O(K^3)$, and the projection, $O(P_zq)$ per target. Because the retained modes are nested, one decomposition to the largest requested $q$ gives every smaller rank by cumulative sums. Recomputing it for each rank, as naive scripts do, throws this away. “PCA” by itself does not specify centring, weighting, reference selection, masking or locality, and these choices change performance more than the eigensolver.

RDI and the JWST era

JWST’s coronagraphs see a much more stable PSF than ground-based adaptive optics, but the telescope can roll by only a small angle in one visit, so reference-star subtraction carries most of the load. JWST’s first exoplanet images, of the known planet HIP 65426 b, used rolls separated by about 10 degrees (about 14 degrees was the largest available in the same epoch), a reference star HIP 68245, and seven filters from about 2 to 16 micrometres. At 1 arcsec the contrast reached about $10^{-5}$ at 4.4 µm and $2\times10^{-4}$ at 11.3 µm (Carter et al. 2023). The paper’s often-quoted sensitivity to about 0.3 Jupiter masses beyond roughly 100 au is a model-based limit on additional companions; the planet itself has an estimated mass of $7.1\pm1.2\,M_{\rm J}$ from luminosity and evolution models.

Two later JWST results show what reference subtraction can and cannot do. The cold, mature giant $\epsilon$ Indi Ab was imaged with MIRI using RDI with eight PCA modes, cross-checked with non-negative matrix factorisation. Its dynamical mass is $6.31^{+0.60}_{-0.56}\,M_{\rm J}$ and its atmosphere fits about 275 K. The Methods section identifies a negative residual caused by a companion of the reference star, a concrete instance of library contamination (Matthews et al. 2024). TWA 7 used nine reference-star dither positions, about an hour on the reference and four hours in total on the target, and reported evidence for a roughly 0.3-Jupiter-mass object at 1.5 arcsec and about 320 K. Its title says “evidence for”, and its mass depends on models (Lagrange et al. 2025).

Even with a planet-free reference library, a target planet can be partly fitted by stellar basis vectors. More references are not automatically better: mismatched modes fit irrelevant variation. Selecting references after looking at a candidate creates selection trials.

Self-subtraction and forward modelling

Two distinct losses afflict every subtraction method. Over-subtraction projects planet light onto stellar modes even when the references are clean. Self-subtraction happens when planet light enters the references and changes the modes themselves. Pueyo’s 2016 KLIP-FM writes the reference set as $R(f)=R_0+fA$ for a planet of flux $f$, expands the eigenimages as $Z(f)=Z_0+f\,\delta Z+O(f^2)$, and obtains the planet response to first order:

The first term is fixed-basis attenuation; the second is the response of the basis. Pueyo derived the eigenvector perturbations analytically and extended them to combined angular and spectral processing (Pueyo 2016). Bright sources, aggressive truncation or nearly degenerate eigenvalues break the linear approximation.

Jean-Baptiste Ruffio and collaborators built on this a forward-model matched filter (FMMF) that uses the processed template at every position. Processing 330 Gemini Planet Imager Exoplanet Survey (GPIES) datasets and comparing methods at a matched empirical false-positive rate, they found a median contrast gain of about 25% over a Gaussian matched filter for T-type spectral templates, sometimes more than a factor of two, and about 50% over simple cross-correlation, equivalent to roughly 2.3 times less exposure time under square-root-of-time scaling (Ruffio et al. 2017). The comparison is unusually valuable because it fixes false positives and completeness, not only residual noise.

The alternative is brute force: the negative fake companion method inserts a negative planet model before the reduction and adjusts its position, flux and spectrum until the residual vanishes. It handles nonlinear response but needs many full reductions inside an optimiser or sampler. Analytical forward models amortise those reductions where their perturbation theory holds.

Statistical detectors: ANDROMEDA, PACO and TRAP

A second family of methods skips the “clean image” stage and fits a statistical model to the data directly.

ANDROMEDA, from maximum-likelihood work by Laurent Mugnier and collaborators and validated on VLT/NaCo data by Faustine Cantalloube and colleagues in 2015, forms differences of frames $i$ and $j$, $\Delta_k=I_i-\gamma_kI_j$, whose planet signature is a positive and a negative copy of the PSF, $p_k=h_i(\boldsymbol\rho)-\gamma_kh_j(\boldsymbol\rho)$. Under diagonal Gaussian noise the flux estimate is the matched filter $\widehat f=\sum p_k\Delta_k/\sigma_k^2\big/\sum p_k^2/\sigma_k^2$. The paper reports about 180 seconds for 612 images of $150\times150$ pixels on a multi-core Xeon workstation, compared with about a second for a PCA residual image but about ten minutes for an injection-based detection limit (Cantalloube et al. 2015). Different products, different hardware: the comparison mainly shows that calibration can cost more than subtraction.

PACO, introduced by Olivier Flasseur, Loïc Denis, Éric Thiébaut and Maud Langlois in 2018, models each small patch of the image as a local mean plus Gaussian noise with a locally learned spatial covariance,

where $S_x$ is the sample covariance and $\rho_x$ a data-driven shrinkage weight that trades estimation noise against lost correlation. Following each trial sky position through the detector as the field rotates, PACO accumulates

the same $A$ and $B$ as before (Flasseur et al. 2018). Its fast variant learns the local statistics once and precomputes $w_x=\widehat C_x^{-1}h_x$, $c_x=h_x^Tw_x$ and $d_{x,t}=w_x^T(r_{x,t}-\widehat m_x)$; detection then becomes sampling the $c$ and $d$ maps along the planet’s track and summing. For 96 frames of $1024^2$ pixels with 49-pixel patches, in MATLAB on 24 Xeon cores, the fast version took about two minutes against three hours for the complete one, a reduction of roughly 90 times obtained by reuse alone, before any GPU. PACO ASDI extended the model to integral-field spectrographs in 2020 (Flasseur et al. 2020). On 24 SPHERE targets Antoine Chomez and colleagues found IRDIS contrast improvements of about seven over TLOCI and five over PCA at 0.5 arcsec, while noting that the competing "$5\sigma$" curves did not share a false-alarm rate (Chomez et al. 2023).

TRAP, by Matthias Samland, Jeroen Bouwman, David Hogg and collaborators (2020 preprint, 2021 paper), turns the regression around. A rotating planet produces a predictable light curve in each pixel it crosses. Other pixels, far enough away not to contain the planet, share the instrumental causes. Their light curves become regressors:

with $U$ the temporal principal components of the reference pixels, $\mathbf s_i$ the planet’s light curve in pixel $i$, and nuisance and planet fitted simultaneously. Because exclusion is spatial rather than temporal, every frame can be used, avoiding ADI’s protection gap. For 51 Eridani the gain reached about six at $2\lambda/D$ relative to ANDROMEDA with $1\lambda/D$ protection, and for $\beta$ Pictoris using the native four-second sampling raised planet S/N by about four relative to the spatial model (Samland et al. 2021). Averaging frames before reduction can irreversibly discard information about changing systematics. TRAP’s cost example, 256 frames within about 550 mas, took about 120 minutes on one laptop core and 40 on four, scaling roughly as $T^{1.5}$.

Surveys and what they found

The GPIES analysis of its first 300 stars detected six planets and three brown dwarfs, including previously known objects; 51 Eridani b and HR 2562 B were the new discoveries in that subsample. It reduced every dataset with 30 KL modes from the 150 most correlated eligible frames and a minimum combined angular-spectral motion of 0.7 pixel, used thresholds of $6\sigma$ outside 0.3 arcsec, and raised them to $8\sigma$ inside for T-type reductions because close-in residuals produced excess false positives. It inferred an occurrence of $9^{+5}_{-4}\%$ for 5–13 Jupiter-mass planets at 10–100 au around its higher-mass hosts (Nielsen et al. 2019). The COBREX reanalysis of 400 GPIES stars with PACO improved minimum detectable masses by up to a factor of two and produced two promising candidates but no confirmed new companion, with an occurrence of $1.7^{+0.9}_{-0.7}\%$ over similar ranges for its sample (Squicciarini et al. 2025). Host samples, selection and evolutionary models differ; the two numbers are not a change in the universe.

SHINE’s first demographic analysis used 150 stars and 13 substellar detections (Vigan et al. 2021). Its final analysis processed 650 datasets of 400 stars uniformly with PACO ASDI and classified more than 3,500 physical sources, most of them background stars, leaving 24 additional promising candidates. Those are not new planets (Chomez et al. 2025).

The NASA Exoplanet Archive credits imaging with 98 confirmed planets through mid-October 2026, against more than 4,700 for transits. The numbers are small because young, self-luminous giants on wide orbits are rare. That is also why imaging’s value lies in what it measures, spectra and orbits, as much as in what it counts.

Stacking along orbits

Why motion matters

Every method above assumes the planet stays put on the sky during an observing sequence, and for young giants tens of au from their stars that is true. It stops being true as searches move inward, to longer integrations, larger telescopes and shorter wavelengths. For a face-on circular orbit of radius $a$ (au) around a star of mass $M_\star$ (solar masses) at distance $d$ (parsecs), the sky motion is

Jared Males, Andrew Skemer and Laird Close expressed it in units of the PSF width per day (Males et al. 2013):

For an Earth-like orbit around a Sun-like star at 10 pc, seen with an 8-metre telescope in the H band (1.6 µm), this gives about 0.04 FWHM per day: about 1.3 PSF widths in a month and some 15 in a year (an illustrative calculation). Nights a month apart cannot simply be coadded. For the NEAR campaign on $\alpha$ Centauri A, Hervé Le Coroller and collaborators give 12.8 mas per day for a face-on circular orbit at 1 au, and a total displacement over the campaign of more than 350 mas, comparable to the mid-infrared PSF (Le Coroller et al. 2022). Larger apertures make the problem worse: better resolution increases the number of distinguishable trajectories even for the same physical orbit.

Males and collaborators noted the connection to moving-object shift-and-stack, and cited Skemer and Close’s use of the known orbit of Sirius B to combine infrared observations as a precursor (a recovery of a known white dwarf, not a discovery). In 2015 Males, Ruslan Belikov and Eduardo Bendek proposed orbital differential imaging, which treats orbital motion as a further kind of diversity, useful for building planet-free references as well as for coaddition (Males et al. 2015).

The orbital signal model

For a single planet in the two-body approximation, the orbit is $\boldsymbol\theta=(a,e,i,\Omega,\omega,\tau)$ with mean motion $n=\sqrt{G(M_\star+M_p)/a^3}$. At epoch $t_j$ solve Kepler’s equation,

form the orbital-plane coordinates $X_j=a(\cos E_j-e)$ and $Y_j=a\sqrt{1-e^2}\sin E_j$, rotate by $R_z(\Omega)R_x(i)R_z(\omega)$, project onto the sky, divide by $d$, and apply each epoch’s plate scale, north angle, distortion and stellar centre to get the predicted detector position $\boldsymbol x_j(\boldsymbol\theta)$ (Nowak et al. 2018; Dallant et al. 2023). The reduced data at that epoch are modelled as

where $\boldsymbol h$ is the response to a unit-flux planet after the whole reduction, negative lobes and all. Ascending-node conventions and east–west handedness are easy places to make a plausible-looking error.

Notice what is absent. A star and its planet share the same parallax and proper motion, so in star-centred coordinates the observer’s motion cancels exactly and the planet’s track is pure Kepler. This is the reverse of the distant-Solar-System case of Chapter 20, where the apparent track of a body at hundreds of au is mostly the reflection of Earth’s own orbit.

K-Stacker

At the October 2015 Haute-Provence meeting “Twenty years of giant exoplanets”, Le Coroller, Mathias Nowak, Luc Arnold, Kjetil Dohlen, Thierry Fusco, Jean-François Sauvage and Arthur Vigan presented Keplerian-Stacker: a blind search over the orbit of a planet invisible in every individual image. Their early demonstration took about ten hours on 190 processors (Le Coroller et al. 2015). Nowak and collaborators’ 2018 paper gave the reproducible design. For a trial orbit, predict one position per epoch, measure the background-subtracted aperture signal $s_j$ there, estimate the annular noise $n_j$ at that radius, and score

The 2018 grid had six orbital dimensions with stellar mass and distance fixed: 10 values of $a$ (0.9–7.5 au for the example star), 10 of $e$ (0–0.8), 100 of the periastron time, 35 each of $\Omega$ and $\omega$, and 7 of $i$. That is 85,750,000 orbits. The 100 best grid points were refined by BFGS and the recombined images inspected. On a cluster with 50 cores, an estimated 150 GFLOP/s in aggregate, each ten-image sequence took 10–15 hours for the grid plus 1–2 hours of refinement (Nowak et al. 2018).

The blind test was honest about its scope. Fifty simulated ten-image sequences were analysed without knowing which contained planets. Of the 35 with injected planets, 25 yielded the correct candidate, nine were missed and one produced a false candidate elsewhere; all 15 empty sequences were correctly declared empty. The simulations contained atmospheric residuals but not realistic ADI processing or long-lived instrumental speckles. In 2020 the team moved injections into real SHINE observations, before PCA/ADI reduction so that self-subtraction was included. Few stars had ten repeated epochs, so sequences were assembled from 50 images of similar stars, and repeated injection gave 30 test runs, six of them empty. Every planet whose ideal recombined S/N exceeded 6.5, about 2 per observation, was recovered. The empty and faint cases produced ambiguous candidates with scores around six, some aligning a real weak source in only a few epochs. The same paper recovered the known planets $\beta$ Pictoris b and HD 95086 b in their own repeated observations (Le Coroller et al. 2020).

calibrate images, stellar centres, PSFs, masks and epoch coordinates
precompute per-epoch flux maps and noise-variance maps
for each orbit theta in the grid:
    signal = variance = 0
    for each epoch j:
        x = project_kepler_orbit(theta, t[j], calibration[j])
        if x is unmasked:
            signal   += interpolate(flux_map[j], x)
            variance += interpolate(noise_map[j], x)
    keep (theta, signal / sqrt(variance)) in a bounded top list
cluster survivors by predicted image track; refine each continuously
calibrate the maximum with null searches and injections

With precomputed maps the cost is about $O[TN(J+k)]$ for $T$ orbits, $N$ epochs, $J$ Kepler iterations and $k$ interpolation samples. No trial ever constructs a shifted image.

The same machine as Chapter 20

Orbit stacking of exoplanets and orbit-space stacking of Solar-System bodies are the same algorithm with different coordinates. Both search for a source invisible in every image by adding evidence along physically allowed trajectories before any detection is required. Both compute, for each image, a pair of matched-filter maps; both sample those maps at predicted positions and form one normalised statistic; both measure template mismatch with the same metric. The parallel is exact enough to be written as a dictionary.

In heliostack, $\xi_k$ is the image divided by its variance and cross-correlated with the PSF, and $\zeta_k$ the inverse variance correlated with the squared PSF; summing them along a trajectory gives $\Psi$ and $\Phi$, with best-fit flux $\Psi/\Phi$ and significance $\nu=\Psi/\sqrt\Phi$ (Napier et al. 2026). PACOME’s $b$ and $a$ maps are the same two quantities, built with PACO’s local covariance instead of a diagonal variance. The mismatch metric derived for multi-year digital tracking, $g=J^T\eta J$ with $J=\partial\boldsymbol x/\partial\boldsymbol\theta$ and $\eta_k=c_kI_2/(4b_k^2)$ (Geringer-Sameth et al. 2025), is the expansion of the overlap formula below. Alex Parker and J. J. Kavelaars’s 2010 construction of near-minimal shift-vector banks from synthetic physical orbits (Parker & Kavelaars 2010) is the right model for exoplanet orbit banks too. The closest Solar-System analogue of all is satellite stacking around a bright primary, such as Burkhart, Ragozzine and Brown’s search for faint moons of Haumea with nonlinear host-relative shifts (Burkhart et al. 2016): host glare and primary-relative orbits make it nearly the exoplanet problem.

The differences are real, and they decide where the difficulty lies. Solar-System stacks have many exposures with nearly independent noise, so the false-alarm problem is mostly one of enormous trial counts. Exoplanet stacks have ten or twenty epochs of correlated, heavy-tailed residuals, so a single surviving speckle can dominate a trajectory’s score. Exoplanet brightness also changes along the orbit: reflected light follows $F_p/F_\star=A_g\Phi(\alpha)(R_p/r)^2$, with geometric albedo $A_g$, phase function $\Phi$ and star–planet distance $r$, and polarised light adds a further phase dependence. And masks, the inner working angle and coronagraph boundaries remove whole epochs from some trajectories. A masked epoch is missing information, not a measurement of zero flux.

Mismatch, trials and what a score means

For a circular Gaussian PSF of standard deviation $s$, two correctly shaped templates displaced by $\delta$ overlap by $\exp[-\delta^2/(4s^2)]$. With normalised statistical weights $w_j$, a trial orbit $\boldsymbol\theta+\Delta\boldsymbol\theta$ retains a fraction

of the matched-filter mean. Expanding to second order gives the metric $1-\eta\simeq\Delta\boldsymbol\theta^Tg\,\Delta\boldsymbol\theta$ with $g=\sum_jw_jJ_j^TJ_j/(4s_j^2)$. It measures distinguishability in the observed image sequence, not distance between orbital elements. Over a short arc many elements barely move any predicted position; a uniform grid in all six or seven elements wastes most of its templates there, while near periastron the same grid can undersample sharply bending tracks. NEAR quantified the stakes: ignoring a seven-pixel displacement cost about 15% of the score, an eleven-pixel one about 25% (Le Coroller et al. 2022).

Trials are the other half. A single standard-normal test exceeds 5 with probability $2.9\times10^{-7}$. If the 85,750,000 orbits of the 2018 grid were independent, about 25 would exceed 5 by chance alone (an illustrative calculation). They are not independent, since neighbouring orbits share nearly all their samples, and the observed false-candidate rate was much lower, partly because people also inspected morphology. Neither fact shows that a universal $5\sigma$ threshold is safe (Nowak et al. 2018). The honest procedure calibrates the maximum statistic of the complete search, including refinement, masks, spectral alternatives and candidate selection, on noise realisations that preserve the real correlations. The same lesson recurs in pulsar and gravitational-wave searches (Chapter 38).

PACOME: likelihood maps and a seven-dimensional search

Jules Dallant, Maud Langlois, Olivier Flasseur and Éric Thiébaut published PACOME in 2023: PACO’s statistical model extended across epochs, with a fully automated orbital search written in Julia. With precision matrix $W$, unit-flux PSF $h$ and frames $k$ within each epoch $j$ and channel $\lambda$, it builds

so that the local flux estimate is $b/a$ with variance $1/a$. Allowing an independent non-negative flux for every epoch and channel and profiling it out gives

twice the log-likelihood improvement. The seventh search dimension is a Kepler constant combining stellar mass and distance, so that dynamical scale uncertainty propagates. A regular grid is followed by refinement of the best seeds with the VMLM-B quasi-Newton method, using analytic derivatives of the orbit and of the spline interpolation (Dallant et al. 2023).

The compute numbers are the most important in this section. For ten epochs, two channels and $10^9$ trial orbits, PACOME took about twenty minutes on twelve cores of an Intel Xeon E5-2620 v3, explicitly without GPUs: about 9 µs per orbit and about 50 µs for a cost-plus-gradient evaluation. The upstream PACO reduction of one 500-frame, $1024\times1024$, two-channel hypercube took about an hour on ten threads. The search maps for ten epochs and two channels occupy about 168 MB in float32. Larger studies tested about $1.7\times10^{11}$ orbits for injected sources (44 hours on twelve threads) and $2.05\times10^{10}$ to fit the four known HR 8799 planets (22 hours). Combining 23 IRDIS and 15 IFS observations of HR 8799, the known planets reached combined scores of 231 to 530.

$\mathcal C$ is not a Gaussian S/N. Under the ideal null of independent standard-normal $b/\sqrt a$ terms, the number of positive terms is binomial and each squared amplitude is independent of its sign, so

a mixture of chi-square distributions with $\chi^2_0$ a point mass at zero. With $M=20$ terms and PACOME’s own per-trial budget of $0.1/N_{\rm orb}=10^{-10}$ for $10^9$ orbits, the threshold is $\sqrt{\mathcal C}\approx8.5$, against 6.4 for a single signed Gaussian statistic at the same tail probability (my calculation from the mixture). Freedom to vary the flux from epoch to epoch has a price in threshold. PACOME itself estimated the distribution numerically and fitted extreme quantiles, with corrections for non-unit variance; correlations, refinement and masking still need empirical tests.

Its failures are as instructive as its successes. In the outer HR 8799 search a high-scoring trajectory drew 41% of its score from two consecutive observing dates, plausibly an instrumental artefact; in the inner search several spectral alternatives crossed the nominal threshold while being dominated by a few epochs. Independent per-epoch fluxes tolerate real variability but also reward a path that strings together isolated positive noise excursions. A common flux suppresses that failure but rejects real sources whose brightness changes with phase. William Thompson and collaborators’ HR 8799 search took the common-flux route, fitting one brightness across 16.75 hours of Keck/NIRC2 data with Hamiltonian Monte Carlo. Its tentative inner candidate, consistent with 4–7 Jupiter masses at 4–5 au, had a marginalised photometric S/N of 6.9 and a log Bayes factor of 18, yet the authors judged its Gaussian-equivalent false-alarm level as only better than about three sigma once non-Gaussian noise and small samples were considered (Thompson et al. 2022). Their Octofitter framework generalises this to images, astrometry, radial velocities and interferometry; in a restricted circular face-on test, a source at per-image S/N 0.5 in 100 monthly epochs was recovered near total S/N five (Thompson et al. 2023).

Campaigns: $\alpha$ Centauri, $\epsilon$ Eridani and the meaning of a candidate

NEAR (New Earths in the Alpha Cen Region), an ESO–Breakthrough Initiatives experiment, used VISIR on the VLT with adaptive optics, a vortex coronagraph and rapid chopping between $\alpha$ Cen A and B. Kevin Wagner and colleagues collected about 100 hours in 2019 and kept 76.9 hours of good data. They reported a candidate, C1, near $\alpha$ Cen A at $0.85\pm0.05$ arcsec with flux $1.2\pm0.4$ mJy and a small-sample S/N of about 3.1, and considered a planet, exozodiacal dust and an instrumental artefact as explanations (Wagner et al. 2021). The 2022 K-Stacker reanalysis used sixteen epochs from late May to late June 2019 at 11.25 µm, masked a persistence strip covering about 36% of the field, fixed the mass and distance at $1.133\,M_\odot$ and 1.34 pc, and searched about $3.5\times10^5$ orbits with $a$ from 0.95 to 2 au. It recovered C1 with an optimised score of about 5.5 and no other candidate. Eighty-four top orbits described essentially the same image track: one object, not 84. For C1 itself orbital compensation gave no significant gain, because its motion over a month was small; for about 15% of the searched orbits it would have raised the score by 20–25%. A synthetic test that spread the same sixteen images over half a year recovered an injected planet at per-image S/N about 1 with score 5.5, which plain coaddition missed. The completeness reached about 30% near eight Earth radii, capped at 60% by the masks (Le Coroller et al. 2022). This was an independent analysis of the same photons, not a confirmation.

In 2025 Charles Beichman, Aniket Sanghi and collaborators imaged $\alpha$ Cen A three times with JWST/MIRI at 15.5 µm. A source S1, about 3.5 mJy and 1.5 arcsec (about 2 au) from the star, appeared in August 2024 at S/N 4–6 depending on processing, but not in February or April 2025. Assuming S1 and C1 are the same object, they drew $10^7$ fitting orbits, kept the roughly 30% that were stable over a million years, and found that 52% of those would have been undetectable in both later epochs. That 52% is a conditional visibility fraction under an assumed identification, prior and brightness, not the probability that S1 is a planet (Beichman et al. 2025; Sanghi et al. 2025). Their public injection-recovery maps show why such reasoning must be done orbit by orbit.

The same median sensitivity at 1.5 arcsec was about 11.5 mJy in August and about 2.7–2.8 mJy in the two later epochs, and each map varies strongly with position angle. A nondetection constrains an orbit only where that orbit puts the planet. The natural extension is a censored likelihood that multiplies, over undetected epochs, the probability $1-c_j(\boldsymbol x_j(\boldsymbol\theta),F)$ of having missed a source of flux $F$ at the predicted position. A coronagraph-hidden position contributes almost nothing; a clean, sensitive one can reject an orbit decisively.

C. Tschudi and collaborators’ 2024 search for $\epsilon$ Eridani b pushed the method to reflected light: 38.5 hours of SPHERE/ZIMPOL polarimetry at 735 nm, twelve nights in four epochs over fourteen months. With stellar mass also free, the grid had $6\times21\times48\times220\times110\times40\times110\approx6.4\times10^{11}$ orbits. A grid score reached 7.4 and refinement typically added 0.4, yet the authors claimed no detection: several competing tracks converged on a strong instrumental speckle ring. The combined polarised-contrast limit near one arcsecond was about $8\times10^{-9}$, about twice as deep as a single epoch (Tschudi et al. 2024).

Learning the speckles

Machine learning in high-contrast imaging is best understood as another way to model the nuisance, and its strongest results so far are improvements in injection-recovery and archival recoveries rather than a population of new planets (Chapter 22 treats ML detection in moving-object surveys).

SODINN, by Carlos Alberto Gomez Gonzalez, Olivier Absil and Marc Van Droogenbroeck (2018), classifies how a local patch changes as PCA subtraction grows stronger. For annular data $M\in\mathbb R^{T\times P_a}$ it forms residuals $R_k=M-MV_kV_k^T$ at 8–20 ranks spanning cumulative explained variance from about 0.5 to 0.99, derotates them and stacks patches across ranks. The stacking axis is subtraction strength, not time: a planet has a characteristic core and self-subtraction pattern that a speckle does not. Two convolutional LSTM layers classify the stack. Positive examples come from injecting the measured PSF into the raw cube; negatives from source-free patches with reversed rotation. Training took about an hour on one P100. At matched false-positive levels near one to two diffraction widths, the true-positive rate improved by about 2 to 10 times over PCA, a gain in recovery rate, not in contrast (Gomez Gonzalez et al. 2018). Its 2026 noise-adaptive successor reprocessed the SHINE F150 sample, reported median 95%-completeness contrasts near $10^{-5}$ at 250 mas and $10^{-6}$ at one arcsecond, and found 30 additional point-source candidates: 17 already reported elsewhere and 13 new. Follow-up colours and epochs left one, near Smethells 20, as the strongest remaining target. The full workflow took 10–30 hours per sequence and band on CPUs, dominated by the injections needed for calibration (Cantero et al. 2026).

Deep PACO (Flasseur, Bodrito and colleagues, 2023) keeps PACO’s whitening and learns only the remaining discrimination, with a U-Net on whitened, derotated data. Frames are temporally shuffled before training injections, which breaks real-source trajectories while keeping local noise. Across eleven SPHERE/IRDIS datasets, detection contrast improved by up to half a magnitude over PACO, a flux factor of 1.58. The paper’s careful validation is its most useful contribution: training once and injecting afterwards, instead of retraining for each injection, cut the cost of sensitivity estimation by about 2,000 and was valid for detection, but overstated photometric accuracy close to the star by up to a factor of five (Flasseur et al. 2023). The output is a pseudo-probability map, which the authors decline to read as a probability that a planet is present.

MODEL&CO (Théo Bodrito and colleagues, 2024) learns across an archive: 220 curated cubes with 20,462 frames, excluding disks, crowded fields and short sequences. It processes features in speckle-aligned coordinates, whitens them with an estimated covariance, combines them in sky coordinates, and calibrates its output on 100,000 augmented realisations of independent calibration observations to estimate the tail of the null distribution. For limited-rotation sequences it reports contrast gains of eight to ten over PACO. Training takes a few hours on V100 or GTX 1080 Ti GPUs, inference a few seconds (Bodrito et al. 2024).

4S (Markus Bonse and colleagues, 2024) is a linear noise model with constrained information paths. The stellar estimate is $x_tB$, where the weights predicting a pixel are masked near that pixel and regularised with the PSF so the model cannot copy a planet into its own prediction. The loss is evaluated after derotation,

so a constant sky source contributes to the derotated mean rather than to the variance being minimised. Optimisation with L-BFGS takes 10–30 minutes on one A100 or H100. It recovered AF Lep b in VLT/NaCo data from 2011, extending that planet’s astrometric baseline by eleven years (Bonse et al. 2024). The gain came from protecting the signal, not from depth.

Kai Hou Yip and colleagues trained a CNN on GAN-generated stellar residuals (Yip et al. 2019), and a 2025 CNN–transformer reached 100% accuracy on a small synthetic test (Lin 2025), which its author calls a proof of concept. Perfect accuracy on a balanced test of 500 samples says nothing about the rare-tail reliability a blind search needs. The central difficulty of ML here is calibration, not classification: real planets are few, injections dominate the positives, adjacent patches and frames of the same star leak information between training and test sets, and observing conditions shift.

Interferometry: measuring the planet in Fourier space

Interferometric methods distinguish planet from star with observables that cancel instrumental errors, before an image is made.

For an unresolved star and companion with flux ratio $r_\lambda$ and separation $\boldsymbol\theta$, the normalised complex visibility on a projected baseline $\mathbf B$ is

Each telescope adds an unknown piston phase, $\arg\widetilde V_{ij}=\arg V_{ij}+\phi_i-\phi_j$, but the phases cancel around a closed triangle: $\Phi_{ijk}=\arg(\widetilde V_{ij}\widetilde V_{jk}\widetilde V_{ki})$. Roger Jennison introduced closure phase in radio interferometry in 1958. Aperture masking brings it inside one telescope: a mask of $N$ holes gives $N(N-1)/2$ baselines and $(N-1)(N-2)/2$ independent closure phases, so a seven-hole mask has 21 baselines and 15 independent closure phases although 35 triangles can be listed. Closure phases respond only to asymmetric structure, and for faint companions very close to the star the response becomes cubic in baseline times separation, so “inside the diffraction limit” does not mean arbitrarily sensitive. Frantz Martinache’s 2010 kernel phase generalised the idea to a full pupil: in the small-aberration regime, phases are $\Phi=R^{-1}A\varphi+\Phi_{\rm object}$, and any matrix $K$ with $KA=0$ removes the instrumental term. His HST/NICMOS demonstration used 156 pupil samples, 366 Fourier samples and 288 kernel relations, and recovered the known stellar companion of GJ 164 at $88.5\pm3.6$ mas (Martinache 2010). JWST’s aperture masking of HIP 65426 at 3.8 µm reached about 7.8 magnitudes at roughly 70 mas but detected nothing, performed one to two magnitudes worse than predicted, and implicated detector charge migration and calibrator mismatch (Ray et al. 2023).

GRAVITY at the VLT Interferometer has detected planetary light directly. On 28 August 2018, with the four 8-metre Unit Telescopes, one fibre tracked fringes on HR 8799 while another collected light at the position of planet e, about 390 mas away. Stellar leakage and the planet’s coherent flux were separated in a joint fit, yielding a $K$-band spectrum at resolution 500 and astrometric uncertainties of 55 and 190 microarcseconds along the axes of the error ellipse (GRAVITY Collaboration 2019). Its fibre accepts only about 60 mas on the sky, so it must be pointed where a planet is predicted to be. In 2020 Nowak and the GRAVITY team pointed it at the position predicted by the radial-velocity orbit of $\beta$ Pictoris c and confirmed it directly, measuring $a=2.72\pm0.02$ au and $M=8.2\pm0.8\,M_{\rm J}$ (Nowak et al. 2020). Sasha Hinkley and colleagues, guided by an unexplained radial-velocity drift and a proper-motion anomaly, discovered HD 206893 c at $a=3.53^{+0.08}_{-0.07}$ au with $M=12.3^{+1.1}_{-1.2}\,M_{\rm J}$ (Hinkley et al. 2023). Thomas Winterhalder and collaborators followed eight Gaia orbital predictions and detected all eight companions, seven previously unknown and five substellar (Winterhalder et al. 2024). Astrometric targeting (Chapter 34) is what makes a 60-mas fibre a discovery instrument.

Nulling interferometry suppresses starlight by destructive interference. Ronald Bracewell proposed it for exoplanets in 1978. For a two-aperture nuller the dark-output transmission is $T_-(\boldsymbol\theta,\lambda)=\sin^2(\pi\mathbf B\cdot\boldsymbol\theta/\lambda)$, but a stellar disk of angular radius $\alpha_\star$ leaks about $\pi^2B^2\alpha_\star^2/(4\lambda^2)$, so longer baselines both sharpen the fringes and resolve the star. Kernel nulling, proposed by Martinache and Michael Ireland in 2018, cancels low-order instrumental terms in the nulled outputs (Martinache & Ireland 2018). The proposed LIFE space mission would use four collectors; its signal-extraction study combines a rotating, wavelength-dependent modulation template into per-channel $a_\lambda$ and $b_\lambda$ and the statistic $Q=\sum_\lambda\max(0,b_\lambda)^2/a_\lambda$, the same structure as PACOME’s free-flux score. For 31 spectral channels the nominal local five-sigma thresholds are about 87 without positivity and 65 with it. An Earth twin at ten parsecs with four two-metre collectors reaches integrated S/N about 9.7 in 200,000 seconds in simulation (Dannert et al. 2022), and the mission study forecasts about 25–45 rocky habitable-zone planets (Quanz et al. 2022). These are forecasts. The engineering status is sobering: the 2026 NOTT chip for the VLTI’s Asgard instrument measured a raw null of $2.4\pm3$ percent at about 138 K (Garreau et al. 2026), while the instrument’s goal is a post-processed contrast of $10^{-5}$.

Landmark results

Limitations and open problems

Identifiability. With too little angular or spectral diversity, planet light resembles the nuisance model and no algorithm can separate them without extra assumptions or extra data. Forward modelling corrects attenuation when information remains; it cannot recover a signal that is degenerate with the star.

Noise models. Covariance estimates are noisy when the number of frames approaches the patch dimension. Original PACO treats frames as independent though speckles persist; TRAP can miss disturbances absent from its reference pixels. A Gaussian objective does not guarantee calibrated extreme tails.

False alarms. The final SHINE analysis found 17 thresholded false positives inside 1 arcsec in IRDIS, eight of them visually recognisable speckles. Controlled reversed-angle tests extrapolated to at least 260 IRDIS and 162 IFS false positives across the survey, and a full reversed-angle reanalysis was deferred because of its compute cost (Chomez et al. 2025). Calibration, not subtraction, is now where compute would most improve reliability.

Planet or disk. A compact dust feature can survive every algorithm. The 2026 SPHERE study of HD 142527 reports a candidate at about 0.09 arcsec (14 au), of 12–50 Jupiter masses if it is a companion, and keeps a disk explanation open (Tran et al. 2026). NEAR’s reanalysis showed that an injected disk edge can attract an orbit search to a point-source-like solution.

Mass is a model output. Converting contrast to mass needs age, distance, atmosphere and cooling models. A factor-of-two gain in minimum detectable mass is not a factor-of-two gain in photon sensitivity.

Short arcs. Many orbits fit a short arc almost equally well. The spread of the top grid solutions depends on the grid, and local Hessian errors miss disconnected modes. Calibration errors in mass, distance, north angle or plate scale can be absorbed into biased orbital elements, which helps detection and hurts interpretation (Le Coroller et al. 2020).

Where the room is

Physical limits are photon noise from leaked starlight, diffraction, coronagraph throughput, and the diversity actually acquired: rotation, bandwidth, epochs. Engineering limits are repeated decompositions, redundant covariance learning, resampling, memory traffic and injection loops. Between them sit model limits, such as discarded temporal structure or false stationarity, that reduce sensitivity even with unlimited FLOPs.

Reuse has already been harvested once. Fast PACO’s 90-fold gain came from learning local statistics once and turning detection into sampling and summing precomputed maps. A GPU port that tiles patch statistics, batches small Cholesky factorisations and fuses interpolation with accumulation could plausibly gain another order of magnitude over the historical MATLAB code, but that is unmeasured, and the 90-fold step must not be claimed again. For KLIP, one decomposition to the largest rank serves all ranks, and frames sharing eligible references can share products. For TRAP, neighbouring trial positions have overlapping exclusion masks, so cached decompositions with low-rank updates could replace a fresh factorisation per position. Each approximation must be validated on faint injected planets, not on the norm of the subtracted starlight. KLIP’s products map well onto dense GEMM and tensor cores; PACO needs many small factorisations and TRAP many small least-squares fits; FP16 throughout is dangerous because bright structure is subtracted to reveal faint residuals and normal equations square the conditioning (Chapter 40).

The orbit search is now a small problem. PACOME’s benchmark, $10^9$ orbits over ten epochs and two channels, implies $2\times10^{10}$ map samples. A bilinear read of two maps is eight float32 values per sample, about 640 GB of logical traffic before cache reuse, and the maps themselves total 168 MB. An H100’s memory bandwidth of about 3.35 TB/s would stream that in roughly 0.2 s, and the $10^{10}$ Kepler solves amount to about $10^{12}$ floating-point operations even at a hundred per solve, a small fraction of a second of GPU arithmetic. Even allowing an order of magnitude for irregular gathers and transcendental functions, the search should take seconds on one GPU rather than twenty minutes on twelve cores (my estimate, unmeasured). The structural trick is the one in the diagram: for fixed $(a,e,\tau,M)$ the eccentric anomalies and orbital-plane coordinates are shared by every orientation, which is just a $2\times3$ projection, so ordering the bank by intrinsic track removes most Kepler solves and keeps map reads local. High eccentricity needs a robust solver; a fixed low iteration count silently loses exactly the orbits being accelerated. The kernel is an irregular gather, not a matrix multiply, and labelling it otherwise does not make it one. The heliostack and LLNL implementations described in Chapter 20 are the natural starting point.

Banks should be built in observable space. The metric $g=J^T\eta J$ says that short arcs need few templates along their poorly constrained directions and long, curved arcs need many. Generating physically allowed orbits and covering their projected tracks, as Parker and Kavelaars did for Solar-System bodies, avoids a rectangular grid in elements that is simultaneously wasteful and too coarse near periastron (Chapter 38). For circular face-on orbits a polar transform turns the search over starting phase into a circular correlation computable by FFT for each radius and frequency (Chapter 39), useful as a seed generator, not as an exact general method. Hierarchical pruning must carry a bound on the fine statistic or a measured completeness loss; rejecting a trajectory because one epoch looks weak removes the very planets orbit stacking was built to find.

Spend the saved compute on the null. If the search takes seconds, thousands of null searches on reversed-angle or shuffled data, and injection campaigns over orbital populations, become affordable. The analytic chi-square mixture gives the ideal null of the free-flux statistic exactly; realistic nulls still need simulation. A hierarchical flux model between “constant” and “free per epoch” could keep tolerance to real variability without rewarding paths through isolated noise peaks, but it must be judged by completeness at a fixed global false-alarm rate, not by its best score. Censored likelihoods that use the published two-dimensional sensitivity maps would make nondetections count properly. External information from radial velocities and Gaia accelerations shrinks the trajectory bank and lowers thresholds, provided it is declared before the search rather than chosen after a candidate appears.

The missing benchmark. No common test measures completeness, false-alarm rate, photometric and astrometric bias, and end-to-end compute across heterogeneous GPI, SPHERE and JWST data. Building one, with frozen search regions, spectral banks and reference policies, is probably the single most useful contribution a methods group could make; Chapter 42 returns to it. Hard limits remain: photons not collected, epochs at which the coronagraph hides the planet, and a planet degenerate with the speckles. No GPU recovers information absent from the data.

References

1. Amara, A., Quanz, S. P. (2012). PynPoint: an image processing package for finding exoplanets. MNRAS 427, 948. arXiv:1207.6637
2. Beichman, C., Sanghi, A., Mawet, D. et al. (2025). Worlds Next Door: A Candidate Giant Planet Imaged in the Habitable Zone of α Cen A. I. Observations, Orbital and Physical Properties, and Exozodi Upper Limits. Accepted manuscript, NASA-hosted. PDF
3. Bodrito, T. et al. (2024). MODEL&CO: exoplanet detection in angular differential imaging by learning across multiple observations. Preprint. arXiv:2409.17178
4. Bonse, M. J. et al. (2024). Use the 4S (Signal-Safe Speckle Subtraction): explainable machine learning reveals the giant exoplanet AF Lep b in high-contrast imaging data from 2011. Preprint. arXiv:2406.01809
5. Burkhart, L. D., Ragozzine, D., Brown, M. E. (2016). A deep search for additional satellites around the dwarf planet Haumea. arXiv preprint. arXiv:1605.03941
6. Cantalloube, F., Mouillet, D., Mugnier, L. M. et al. (2015). Direct exoplanet detection and characterization using the ANDROMEDA method: performance on VLT/NaCo data. A&A 582, A89. arXiv:1508.06406
7. Cantero, C., Sabalbal, M., Absil, O. et al. (2026). New substellar candidates identified through deep learning in the F150 sample of the large-scale SHINE direct imaging survey. Preprint. arXiv:2605.23700
8. Carter, A. L., Hinkley, S., Kammerer, J. et al. (2023). The JWST Early Release Science Program for Direct Observations of Exoplanetary Systems I: high contrast imaging of the exoplanet HIP 65426 b from 2–16 µm. ApJL. arXiv:2208.14990
9. Chomez, A., Lagrange, A.-M., Delorme, P. et al. (2023). Preparing an unsupervised massive analysis of SPHERE high contrast data with the PACO algorithm. A&A 675, A205. arXiv:2305.08766
10. Chomez, A., Delorme, P., Lagrange, A.-M. et al. (2025). The SPHERE infrared survey for exoplanets (SHINE): complete observations, data reduction and analysis, detection performances, and final results. A&A. arXiv:2501.12002
11. Dallant, J., Langlois, M., Flasseur, O., Thiébaut, É. (2023). PACOME: optimal multi-epoch combination of direct imaging observations for joint exoplanet detection and orbit estimation. A&A 679, A38. arXiv:2309.08679
12. Dannert, F. A. et al. (2022). Large Interferometer For Exoplanets (LIFE): II. Signal simulation, signal extraction and fundamental exoplanet parameters from single epoch observations. A&A. arXiv:2203.00471
13. Flasseur, O., Denis, L., Thiébaut, É., Langlois, M. (2018). Exoplanet detection in angular differential imaging by statistical learning of the nonstationary patch covariances. A&A 618, A138. HAL
14. Flasseur, O., Denis, L., Thiébaut, É., Langlois, M. (2020). PACO ASDI: an algorithm for exoplanet detection and characterization in direct imaging with integral field spectrographs. A&A 637, A9. HAL
15. Flasseur, O. et al. (2023). deep PACO: combining statistical models with deep learning for exoplanet detection and characterization in direct imaging at high contrast. Preprint. arXiv:2303.02461
16. Garreau, G. et al. (2026). Asgard/NOTT: cryogenic characterization of the mid-infrared chip. Preprint. arXiv:2607.20191
17. Geringer-Sameth, A., Golovich, N., Iwabuchi, K. (2025). Multi-year stacking searches for solar system bodies. arXiv preprint. arXiv:2509.25428
18. Gomez Gonzalez, C. A., Absil, O., Van Droogenbroeck, M. (2018). Supervised detection of exoplanets in high-contrast imaging sequences. A&A 613, A71. arXiv:1712.02841
19. Gratton, R., Zurlo, A., Le Coroller, H. et al. (2020). Searching for the near infrared counterpart of Proxima c using multi-epoch high contrast SPHERE data at VLT. A&A 638, A120. arXiv:2004.06685
20. GRAVITY Collaboration, Lacour, S., Nowak, M. et al. (2019). First direct detection of an exoplanet by optical interferometry: astrometry and K-band spectroscopy of HR 8799 e. A&A. arXiv:1903.11903
21. Hinkley, S. et al. (2023). Direct discovery of the inner exoplanet in the HD 206893 system. A&A. arXiv:2208.04867
22. Lafrenière, D., Marois, C., Doyon, R., Nadeau, D. et al. (2007). A new algorithm for point spread function subtraction in high-contrast imaging: a demonstration with angular differential imaging. ApJ 660, 770. arXiv:astro-ph/0702697
23. Lagrange, A.-M., Wilkinson, C., Mâlin, M. et al. (2025). Evidence for a sub-Jovian planet in the young TWA 7 disk. Nature. doi:10.1038/s41586-025-09150-4
24. Le Coroller, H., Nowak, M., Arnold, L. et al. (2015). K-Stacker, a new way of detecting and characterizing exoplanets with high contrast imaging instruments. Conference proceedings. arXiv:1510.06331
25. Le Coroller, H., Nowak, M., Delorme, P. et al. (2020). K-Stacker, an algorithm to hack the orbital parameters of planets hidden in high-contrast imaging: first applications to VLT SPHERE multi-epoch observations. A&A 639, A113. arXiv:2004.12878
26. Le Coroller, H., Nowak, M., Wagner, K. et al. (2022). Efficiently combining α Cen A multi-epoch high-contrast imaging data: application of K-Stacker to the 80 hrs NEAR campaign. A&A 667, A142. arXiv:2207.02137
27. Lin, Y.-C. (2025). Machine learning revolution for exoplanet direct imaging detection: transformer architectures. SPIE Proceedings. arXiv:2508.14508
28. Lyot, B. (1939). A study of the solar corona and prominences without eclipses. MNRAS 99, 580. ADS
29. Males, J. R., Skemer, A. J., Close, L. M. (2013). Direct imaging in the habitable zone and the problem of orbital motion. arXiv preprint. arXiv:1304.5853
30. Males, J. R., Belikov, R., Bendek, E. (2015). Orbital differential imaging: a new high-contrast post-processing technique for direct imaging of exoplanets. arXiv preprint. arXiv:1510.02478
31. Marois, C., Lafrenière, D., Doyon, R., Macintosh, B., Nadeau, D. (2006). Angular differential imaging: a powerful high-contrast imaging technique. ApJ 641, 556. arXiv:astro-ph/0512335
32. Marois, C., Macintosh, B., Barman, T., Zuckerman, B. et al. (2008). Direct imaging of multiple planets orbiting the star HR 8799. Science 322, 1348. arXiv:0811.2606
33. Martinache, F. (2010). Kernel phase in Fizeau interferometry. ApJ 724, 464. arXiv:1009.3933
34. Martinache, F., Ireland, M. J. (2018). Kernel-nulling for a robust direct interferometric detection of extrasolar planets. A&A. arXiv:1802.06252
35. Matthews, E. C., Carter, A. L., Pathak, P. et al. (2024). A temperate super-Jupiter imaged with JWST in the mid-infrared. Nature. doi:10.1038/s41586-024-07837-8
36. Mawet, D., Milli, J., Wahhaj, Z. et al. (2014). Fundamental limitations of high contrast imaging set by small sample statistics. ApJ 792, 97. arXiv:1407.2247
37. Napier, K. J., Holman, M. J., Lin, H.-W., Gerdes, D. W., Ruch, T. R. (2026). heliostack: a novel approach to minor planet discovery. Accepted manuscript. arXiv:2604.16731
38. NASA Exoplanet Archive (2026). Planetary Systems table, discovery method counts, retrieved 16 October 2026. NASA Exoplanet Science Institute. exoplanetarchive.ipac.caltech.edu
39. Nielsen, E. L., De Rosa, R. J., Macintosh, B. et al. (2019). The Gemini Planet Imager Exoplanet Survey: giant planet and brown dwarf demographics from 10–100 AU. AJ 158, 13. arXiv:1904.05358
40. Nowak, M., Le Coroller, H., Arnold, L. et al. (2018). K-Stacker: Keplerian image recombination for the direct detection of exoplanets. A&A 615, A144. arXiv:1804.02192
41. Nowak, M. et al. (2020). Direct confirmation of the radial-velocity planet β Pictoris c. A&A. arXiv:2010.04442
42. Parker, A. H., Kavelaars, J. J. (2010). Pencil-beam surveys for trans-Neptunian objects: novel methods for optimization and characterization. PASP. arXiv:1004.3287
43. Pueyo, L. (2016). Detection and characterization of exoplanets using projections on Karhunen–Loève eigenimages: forward modeling. ApJ 824, 117. arXiv:1604.06097
44. Quanz, S. P. et al. (2022). Large Interferometer For Exoplanets (LIFE): I. Improved exoplanet detection yield estimates for a large mid-infrared space-interferometer mission. A&A. arXiv:2101.07500
45. Ray, S. et al. (2023). The JWST Early Release Science Program for Direct Observations of Exoplanetary Systems III: aperture masking interferometric observations of the star HIP 65426 at 3.8 µm. Preprint. arXiv:2310.11508
46. Ruffio, J.-B., Macintosh, B., Wang, J. J., Pueyo, L. et al. (2017). Improving and assessing planet sensitivity of the GPI Exoplanet Survey with a forward model matched filter. ApJ 842, 14. arXiv:1705.05477
47. Samland, M., Bouwman, J., Hogg, D. W., Brandner, W. et al. (2021). TRAP: a temporal systematics model for improved direct detection of exoplanets at small angular separations. A&A 646, A24. arXiv:2011.12311
48. Sanghi, A., Beichman, C., Mawet, D. et al. (2025). Worlds Next Door II: binary star modeling, planet and exozodi search, and sensitivity analysis. Accepted manuscript, NASA-hosted. PDF; data doi:10.5281/zenodo.15748931
49. Soummer, R., Pueyo, L., Larkin, J. (2012). Detection and characterization of exoplanets and disks using projections on Karhunen–Loève eigenimages. ApJL 755, L28. arXiv:1207.4197
50. Squicciarini, V., Mazoyer, J., Lagrange, A.-M. et al. (2025). The COBREX archival survey: improved constraints on the occurrence rate of wide-orbit substellar companions. I. A uniform re-analysis of 400 stars from the GPIES survey. A&A 693, A54. arXiv:2411.06157
51. Thompson, W., Marois, C., Do Ó, C. R. et al. (2022). Deep orbital search for additional planets in the HR 8799 system. arXiv preprint. arXiv:2210.14213
52. Thompson, W., Lawrence, J., Blakely, D. et al. (2023). Octofitter: fast, flexible, and accurate orbit modeling to detect exoplanets. AJ 166, 164. arXiv:2402.01971
53. Tran, T. M. H., Langlois, M., Flasseur, O. et al. (2026). Exploration of the inner region of the system HD 142527. A&A. arXiv:2605.04869
54. Tschudi, C. et al. (2024). SPHERE RefPlanets: search for ε Eridani b and warm dust. arXiv preprint. arXiv:2404.19504
55. Vigan, A., Fontanive, C., Meyer, M., Biller, B. et al. (2021). The SPHERE infrared survey for exoplanets (SHINE). III. The demographics of young giant exoplanets below 300 au with SPHERE. A&A 651, A72. HAL
56. Wagner, K., Boehle, A., Pathak, P. et al. (2021). Imaging low-mass planets within the habitable zone of α Centauri. Nature Communications 12, 922. arXiv:2102.05159
57. Winterhalder, T. O. et al. (2024). Combining Gaia and GRAVITY: characterising five new directly detected substellar companions. arXiv preprint. arXiv:2403.13055
58. Yip, K. H. et al. (2019). Pushing the limits of exoplanet discovery via direct imaging with deep learning. Preprint. arXiv:1904.06155