Direct Exoplanet Imaging Below the Telescope Stability Limit

FLIP is a model-based point-spread-function (PSF) subtraction framework designed for high-contrast observations with the James Webb Space Telescope (JWST) NIRCam coronagraphs. It uses a physical optical model optimized in two stages: first on a reference star to calibrate the model, and then on the science observation to account for changes in pointing and wavefront error. Across four datasets, the method is compared with the standard KLIP algorithm and generally achieves higher signal-to-noise ratios, particularly at small angular separations where telescope instability and imperfect starlight subtraction limit sensitivity.

Model design

FLIP models the stellar PSF by optimizing many free parameters in a physical model that describes both the observed star and the instrument. The models are based on the nominal STPSF optical models used to simulate JWST PSFs for different instruments and filters. Extended Data Fig. 1 presents a schematic equivalent to the physical models used to simulate observations with the JWST NIRCam coronagraphs.

Experiments show that four groups of free parameters are sufficient to reproduce the PSF accurately:

the star’s flux and precise position behind the coronagraph’s focal-plane mask;
the optical path difference (OPD), representing a phase aberration, before the focal-plane mask;
the OPD after the focal-plane mask; and
the response of each detector pixel.

These parameters represent several physical imperfections: telescope pointing errors, imperfect alignment of JWST’s primary mirror segments relative to the closest-in-time OPD measurement, imperfectly known instrument-internal OPD in the NIRCam coronagraphs, and imperfect detector calibration. Other mechanisms can also affect the PSF, including pupil shear, manufacturing defects in the optics and more complex detector physics. However, the authors find that the selected parameters can generally reproduce the effects of these additional systematics and produce high-quality PSF subtractions. A summary of the model parameters is provided in Extended Data Table 1.

Implementation

Following STPSF, FLIP propagates light using Fourier optics in the Fraunhofer approximation. For each observation, STPSF is used to calculate the nominal instrument state and generate a simulated PSF, incorporating the instrument and filter, the OPD measurement closest in time to the reference-star observation, and the host star’s spectrum. This instrument state is then replicated in PyTorch²⁸. Details of the STPSF model construction are given in ref. 27, while implementation details are provided in the Supplementary Methods.

Optimization design

Earlier work using simulated monochromatic data²⁵ found that a nearly exact physical model could be optimized directly on science observations to achieve excellent PSF subtraction without overfitting planetary signals. Tests with real data, however, show that accidental subtraction of planet signals can occur. This is probably because STPSF models are substantially different from real JWST instruments and because FLIP has greater flexibility, with more free parameters making the optimization problem less straightforward.

FLIP therefore uses a two-stage optimization:

Reference-star calibration. The full flexibility of the model is used to optimize all free parameters on a reference star, which is unlikely by design to contain contaminants such as planets or binary companions. This helps compensate for inaccuracies in the optical model and moves the parameters closer to a suitable solution.
Science-image optimization. The calibrated model is then optimized directly on the science observation, but only a reduced set of parameters is varied. This accounts for differences in pointing and OPD between the reference and science observations caused by telescope instabilities.

In both stages, parameters are optimized in separate sub-stages to reduce degeneracies. Stage 1 has three sub-stages:

solve for the star’s position and flux;
solve for both OPDs; and
solve for both OPDs and the pixel responses.

Stage 2 has two sub-stages:

solve for the star’s position and flux; and
solve for the OPD after the focal-plane mask.

Although the OPD before the coronagraph mask is also expected to change as a result of JWST mirror motion, optimization is currently more stable when the OPD after the mask is varied instead. This OPD does not need to be propagated through the coronagraph mask. The stages and their iteration counts are summarized in Extended Data Fig. 2. On a single NVIDIA A30 GPU, stage 1 takes 10–15 minutes and stage 2 takes 10–30 minutes.

Parameter tests and model flexibility

“Leave-one-out” tests illustrate the sources of the contrast improvement (Extended Data Fig. 3). The star’s position has little effect on contrast, although it remains physically important and is included in the model. Its effect may be absorbed by other OPD terms that act as tip-tilt phase aberrations.

OPDs have a strong effect on the PSF because they are directly related to it through a Fourier transform. Detector-response parameters are needed to remove pixel-level errors that cannot be represented by OPDs. This helps explain the improvement over empirical reference-image methods such as KLIP, which cannot account for arbitrary differences between reference and science data during the science-image optimization stages.

Model calibration is also crucial. Removing the calibration stage substantially increases the risk of over-subtracting planetary signals. The tests further show that one OPD plane does not provide sufficient flexibility: in this formulation, the interaction between the two OPD planes separated by the focal-plane mask is important. This distinguishes FLIP from earlier autodiff-based work, particularly ref. 25, which used a less flexible model without model calibration.

Optimizer and learning rates

FLIP uses the Adam optimizer²⁹ with a cosine-annealing learning-rate scheduler and an L1 loss function. Learning rates are selected manually by increasing them across orders of magnitude until the optimization is stable and converges in fewer than approximately 15,000 iterations. Because the model parameters have different numerical scales, each family of parameters is assigned its own learning rate. A coarse search is sufficient to identify the appropriate scale for each parameter.

Extended Data Fig. 4 shows that order-unity changes to the learning rates do not substantially affect the results. Changes by orders of magnitude, however, reveal the scale required for stable and relatively rapid convergence. The approximate learning-rate scales for each parameter are listed in Extended Data Table 1.

Across observations, multiplying all learning rates by an order-unity constant can improve stability and convergence speed. This is probably because the loss-function magnitude depends on stellar brightness. A global scaling factor of four works well for the experiments presented here. Preliminary tests on fainter stars produced better results with scaling factors of one or 0.5.

During the final optimization sub-stage, the global scaling factor is reduced to approximately 1% of its original value to further reduce the risk of overfitting. Results are not strongly sensitive to the precise value, although it can be adjusted for different observations according to the similarity between the reference star and the science target.

Choosing when to stop

To determine when to stop the optimization, the researchers inject artificial planets into the science observation and track their signal-to-noise ratio (SNR) as a function of iteration number. For sensitivity measurements, planets are injected at five to seven separations, with one to three planets per separation. Each artificial planet is generated by simulating an off-axis PSF with the nominal STPSF model at a selected field position and adding it to the science image.

The injected planets have brightnesses between four and 20 times the 5σ fundamental-noise contrast limit, with the value adjusted when necessary according to the observed signal loss. Each planet is injected separately for each FLIP run, and its SNR is measured approximately every 300 iterations.

The SNR generally rises for several iterations and then falls as the algorithm begins to overfit. The optimization is therefore stopped near the iteration at which the injected planets reach an average plateau or peak in SNR. This behavior is consistent for planets injected at similar separations and is illustrated in Extended Data Fig. 5.

For simplicity, the experiments use a single artificial-planet injection to select the stopping point. As with KLIP, the optimal injection strategy depends on the scientific objective. The stopping iteration could instead be chosen using the maximum average SNR across all injections when the goal is to maximize sensitivity across the field, or using one or more injections at a particular brightness and location when optimizing the SNR of a specific planet. Extended Data Fig. 5 shows that PSF-subtraction performance is not strongly dependent on the precise injection choices.

Experiments

The researchers tested FLIP on four datasets: three containing known companions—HIP 65426 in the F444W filter and 51 Eridani in the F410M and F430M filters—and one companion-free target, TYC 5899, from a planet-search survey and observed in F444W.

HIP 65426 and TYC 5899 were observed with the MASK335R coronagraph, while 51 Eridani was observed with the LWBAR coronagraph in the NARROW position offset. The datasets were selected to represent a broad range of JWST NIRCam coronagraphic observations. Together, they include:

three physical systems: HIP 65426, 51 Eridani and TYC 5899;
three filters: F444W, F410M and F430M;
two coronagraphs: the NIRCam 335R round coronagraph and the NIRCam LWBAR bar coronagraph;
real planets, HIP 65416 b and 51 Eridani b, as well as simulated planets in TYC 5899 data;
observations with a dedicated reference star using a 9-POINT-CIRCLE dither pattern for HIP 65426, a dedicated reference star using a 5-POINT-BAR dither pattern for 51 Eridani, and no dedicated reference star for TYC 5899, which instead uses a PSF library from a self-referenced survey;
reference stars observed close in time to the science targets for HIP 65426 and 51 Eridani, and far from the science observations for TYC 5899; and
host stars with spectral types A2V, F0IV and M3 for HIP 65426, 51 Eridani and TYC 5899, respectively.

Standard KLIP PSF subtraction and FLIP were run on all four datasets. Performance was compared using the SNR of recovered companions and detection sensitivity at multiple star–planet separations. KLIP reductions followed the hyperparameters identified as optimal in the relevant publications. In particular, they used reference differential imaging plus angular differential imaging (RDI + ADI; see, for example, ref. 31) on both roll-angle observations of each science target, followed by de-rotation and mean combination. The resulting KLIP subtractions were comparable to those published previously.

For simplicity, FLIP was applied to only one roll angle per dataset. Future work could compare FLIP applied to one roll, FLIP applied to both rolls followed by de-rotation and mean combination, and FLIP applied to both rolls with ADI performed on the resulting images before combining the FLIP + ADI results.

TYC 5899

The original publication on TYC 5899³² presented two epochs: one with a single roll angle because the second roll failed, and another with two roll angles observed four months after the reference-star sequences. This analysis uses the latter epoch for both KLIP and FLIP to demonstrate the effect of substantial drift between the reference and science observations.

Although ref. 32 used all three observations to calculate its contrast curves, this analysis uses only the final two and still obtains comparable contrast performance. At 2.5 arcsec, the contrast is approximately 3 × 10⁻⁶ for KLIP and FLIP. At 1.0 arcsec, it is approximately 4 × 10⁻⁵ for both methods. At 0.5 arcsec, it is approximately 3 × 10⁻⁴ for KLIP and approximately 5 × 10⁻⁴ for FLIP.

Because TYC 5899 has no known companions, planets with a range of separations and brightnesses were injected to compare the two methods under controlled conditions. The FLIP optimization hyperparameters used for the reductions are listed in Extended Data Table 2.

Data preparation

Observations of HIP 65426, 51 Eridani and TYC 5899, together with their reference stars, were processed with JWST pipeline version 1.18.1 and the community-developed spaceKLIP tool, version 2.2.1. The processing steps—including ramp fitting, conversion from detector counts to physical flux units and bad-pixel cleaning—and the observation setups are described in the publications associated with the datasets.^26,31,32

The JWST and spaceKLIP pipelines flag and clean bad pixels, but imperfect cleaning can create more problems for model-based PSF subtraction than for methods such as KLIP. Hot or dead pixels and pixel-response variations cannot be cancelled by subtracting observations that contain the same pixel-level errors. Pixels that remained problematic after automated cleaning were therefore identified by eye and masked during optimization.

Choice of reference stars

HIP 65426 and 51 Eridani were observed with dedicated reference stars: bright stars thought unlikely to host planets, with spectral types similar to those of the science targets, observed shortly before or after the science observations. The FLIP calibration stars were the same reference stars used in the original observations: HIP 68245 for HIP 65426 and HD 30562 for 51 Eridani.

For simplicity, only one dither was used for FLIP model calibration. All dithers were used to construct the principal components for KLIP. Further details are provided in refs. 26 and 31.

TYC 5899 was observed as part of a self-referenced programme, in which only science observations are collected and a large set of observations is combined into a reference library for other science observations. For KLIP, the full PSF library used in ref. 32 was retained to reproduce the published result. For a direct comparison with KLIP, one observation from that library was selected as the FLIP calibration star. The selected observation was HIP 17695, because it has the spectral type most similar to TYC 5899 among the library targets.

PSF-subtraction performance metrics

Contrast curves for KLIP and FLIP were calculated by dividing each image into small concentric annuli and measuring the standard deviation in each annulus. The calculation accounts for small-sample statistics⁵⁰, algorithmic signal loss from overfitting and the throughput of the coronagraph.

Near the measured planets, in the same order as in Fig. 2, algorithmic signal loss for KLIP is approximately 5–10%, 20–70%, 50–80% and 10–40%. For FLIP, it is approximately 5–20%, 20–40%, 30–50% and 40–60%.

All datasets were smoothed with a Gaussian kernel with σ = 1 pixel before standard deviations were calculated, because FLIP residuals are dominated primarily by pixel-to-pixel noise. For HIP 65426, the FLIP contrast curve excludes the central 0.6 arcsec because a bright lobe of the planet’s PSF prevents sampling that region. KLIP uses the KLIP-FM algorithm to subtract the planet before calculating sensitivity. Incorporating planet modelling and subtraction into FLIP could remove this limitation in future work. Details of contrast-curve calculation and calibration are given in refs. 26, 31 and 32.

Fundamental sensitivity

The fundamental sensitivity curve was derived from the ERR extensions of the JWST data products, which track the combined uncertainty in each pixel from read noise, Poisson noise and flat-field noise. A random Gaussian realization was generated for each pixel with a mean of zero and a standard deviation equal to the total noise level. This simulates an image after perfect PSF subtraction.

A contrast curve was then calculated from this noise realization using the same procedure applied to the KLIP and FLIP images, but assuming no planet signal loss. The calculation assumes that the only fundamental noise source is the science observation. This is essentially case 0 from ref. 4. KLIP also receives additional noise from its reference observations, whereas FLIP is model-based and adds no extra fundamental noise. The fundamental-noise curves therefore represent a lower limit, or the most optimistic possible sensitivity.

Because FLIP was applied to only one roll angle per dataset, its noise limit corresponds to the fundamental noise limit of that roll.

SNR maps

To supplement Fig. 2, SNR maps were generated for each PSF-subtracted image (Extended Data Fig. 6). A non-coronagraphic PSF model was cropped to the central 5 × 5 pixels and cross-correlated with the image.

The SNR of each pixel was then calculated within an annulus approximately 0.2 arcsec wide in two steps. First, the noise in each annulus was estimated to identify pixels with SNR > 5. A refined noise estimate was then calculated after excluding those pixels. The final SNR map is the cross-correlation value at each pixel divided by the refined noise estimate.

Stability-limited sensitivity

Figure 3 presents a fiducial stability-limited sensitivity derived by applying KLIP to simulated data. The simulations isolate physical factors that limit KLIP performance. For each observation sequence, 50 realizations were generated with stochastic variations in star–coronagraph misalignment and OPD drift consistent with in-flight JWST characterization.^5,33 No other noise sources, such as photon noise or flat-field noise, were included.

For each science target, STPSF models were generated using the star’s spectral type, and the flux within a central 20-pixel radius was scaled to match the real data. Two rolls separated by 10° in position angle were simulated.

For HIP 65426, nine reference PSFs were simulated using the 9-POINT-CIRCLE dither strategy. For 51 Eridani, five reference PSFs were simulated using the 5-POINT-BAR strategy. For TYC 5899, the same number of reference stars used in the PSF library of ref. 32 was simulated. The relevant dither strategy is described in the JWST documentation.

All reference stars were assumed to have the same spectral type and brightness as the science target. Photon, detector and flat-field noise were omitted. Under these assumptions, the KLIP stability limit can be deeper than the fundamental-noise limit at large separations, as shown in Fig. 3. This is acceptable because the analysis focuses on small separations, where PSF-subtraction inaccuracies rather than fundamental noise generally determine sensitivity.

For each observation, the target-acquisition error was simulated by drawing a random star–coronagraph misalignment from N(μ = 0, σ = 10) mas, consistent with the in-flight target-acquisition accuracy reported in ref. 33. PSFs were generated using the OPD measurement closest in time to each observation, and no significant OPD drift was assumed between reference-star dithers.

For observations requiring a telescope slew—for example, movement from a reference star to the science target—an OPD drift term was added. Its shape was derived by subtracting the two OPD measurements closest in time to the simulated observation. The root-mean-square amplitude of the drift was drawn from N(μ = 4, σ = 1) nm. This is smaller than the approximately 7–9 nm r.m.s. drifts measured over day-long timescales and larger than the approximately 1–2 nm r.m.s. drifts measured over minutes.⁵

For HIP 65426, the bright companion affects PSF subtraction by skewing KLIP, including during mean subtraction and the calculation of Karhunen–Loève modes. The simulations reproduce this effect by injecting a planet with the same contrast, separation and position angle as the real companion.

KLIP was run with ADI + RDI for all targets, using all KL modes for HIP 65426 and 51 Eridani and six KL modes for TYC 5899, following the respective publications. Calibrated stability-limited sensitivity was measured with injection-recovery tests. Planets were injected at three to seven separations and two to three position angles across the 80 × 80 pixel field of view. Their flux was four times the uncalibrated KLIP sensitivity at the injection separation. The central 68% range of the resulting contrast curves is shown as shaded regions in Fig. 3.

TYC 5899 synthetic injections

To compare FLIP and KLIP on TYC 5899, artificial planets with known fluxes were injected across the image and their recovery after PSF subtraction was assessed. Extended Data Figs. 7 and 8 show additional injections that complement Fig. 2.

At each separation, a planet with a flux equal to 20 times the 5σ fundamental-noise limit was injected. All observations were processed identically: FLIP used the same learning rates and iteration counts, while KLIP used the same numbers of annuli and KL modes.

This is a standard way to quantify KLIP sensitivity. Once the annuli and KL modes have been selected—for example, to recover a signal in Fig. 2—artificial planets are injected across the field while those settings remain fixed. The resulting measurements estimate sensitivity for that set of hyperparameters. FLIP uses an analogous procedure, keeping its learning rates and iteration count fixed while injecting planets across the field.

FLIP outperforms standard KLIP in most cases, especially at smaller separations, and produces higher SNR across a diverse range of possible planets. In this dataset, however, FLIP sometimes produces greater signal loss at larger separations. This accounts for the deeper KLIP contrast curve than the FLIP curve at those separations in Fig. 3.

Future work could investigate the difference in signal loss between the methods. At larger separations, the observations are often background-limited rather than contrast-limited, meaning that sensitivity is determined mainly by exposure time rather than by stellar-light subtraction. In the regime most affected by telescope stability and imperfect PSF subtraction—smaller angular separations, which also correspond to more Solar-System-like physical separations—FLIP consistently outperforms KLIP.

Source: www.nature.com