Imaging exoplanets below the telescope stability limit

The FLIP framework

Model design

FLIP models the stellar PSF by optimizing many of the free parameters in a physical model that influence the PSF, including the star being observed and the instrument used for the observation. We base our models on the nominal STPSF optical models, which are used to simulate PSFs in the JWST for a variety of instruments and filters. Extended Data Fig. 1 shows a schematic diagram equivalent to the physical models used to simulate observations with the JWST NIRCam coronagraphs.

Given these models, we find that four groups of free parameters are sufficient for PSF fidelity in our experiments:

• The flux of the star and its precise position behind the focal plane mask of the coronagraph.

• The OPD (a phase aberration) before the focal plane mask of the coronagraph.

• The OPD after the focal plane mask of the coronagraph.

• The response per pixel in the detector.

Physically, these parameters capture the imperfect alignment of the star with the coronagraph (telescope pointing errors), the imperfect alignment of the primary mirrors in the JWST (deformations or tilts relative to the closest-in-time OPD measurement), the imperfectly known OPD internal to specific instruments in the JWST (in our case, the coronagraphs in the NIRCam instrument) and imperfectly known detector calibrations. Although other physical mechanisms are capable of affecting the PSF (such as pupil shear, manufacturing errors in the optics or more complicated detector physics), we find that our choice of free parameters is generally capable of effectively reproducing the effects of these other systematics and produce high-quality PSF subtractions. A summary of our model parameters is shown in Extended Data Table 1.

Implementation details. Following STPSF, our model implements light propagation in the Fraunhofer approximation through Fourier optics. When modelling any specific observation, we use STPSF to calculate the nominal instrument state and produce a simulated PSF, taking into account the instrument and filter, the closest-in-time OPD measurement to the reference star observation and the spectrum of the host star. We take this instrument state and replicate it in PyTorch²⁸. See ref. ²⁷ for details on the STPSF model construction and Supplementary Methods in the Supplementary Information for details on our implementation.

Optimization design

Previous work on simulated monochromatic data²⁵ found that, given a near perfectly accurate physical model, it is possible to optimize it directly on science observations and produce excellent PSF subtraction without overfitting planet signals. However, in our experiments with real data, we find that accidentally subtracting planet signals is indeed possible. This is probably because of STPSF models having substantial inaccuracies relative to real JWST instruments and to our model being much more flexible with more free parameters, thus making the optimization problem less straightforward.

To address this, we implement a two-stage optimization. In stage 1, we use the full flexibility of our model and calibrate it on a reference star (because, by design, it is unlikely to host any contaminants such as planets or binaries), optimizing all of the free parameters described in the ‘Model design’ section. This helps overcome inaccuracies in the optical model and sets the parameters closer to a good solution. In stage 2, we use our calibrated model and optimize only a reduced number of free parameters directly on the science observation. This accounts for different pointing errors and different OPDs between reference and science observations owing to telescope instabilities.

In both stages, we optimize different parameters in separate sub-stages. This is done to reduce degeneracies between different parameters. In stage 1, we have three sub-stages: (1) solving for star position and flux; (2) solving for both OPDs; and (3) solving for both OPDs and pixel responses. In stage 2, we have two sub-stages: (1) solving for star position and flux and (2) solving for the OPD after the focal plane mask. Although the OPD before the coronagraph mask is also expected to change owing to the JWST mirrors, we at present find the optimization to be more stable if we solve for the OPD after the mask instead (as this OPD does not need to be propagated through the coronagraph mask). These stages are summarized and illustrated in Extended Data Fig. 2, with the number of iterations dedicated for each one. On a single NVIDIA A30 GPU, stage 1 takes 10–15 min and stage 2 takes 10–30 min.

We carry out ‘leave one out’ tests to illustrate intuitively the reasons behind our contrast improvement (Extended Data Fig. 3). In the top two rows of Extended Data Fig. 3, we see that the star position has little effect on the contrast: although it is physically important (hence our inclusion), the effect is probably absorbed by other OPD terms in the model as tip-tilt phase aberrations. It is clear that OPDs strongly affect the PSF, as they are related to it directly by means of a Fourier transform. Detector responses are required to remove pixel-level errors that are not the result of and cannot be fit by OPDs. This illustrates the source of the improvements relative to approaches such as KLIP, as methods based on empirical reference images cannot account for arbitrary differences between reference and science data (stages 2.1 and 2.2). The bottom row of Extended Data Fig. 3 shows that model calibration is crucial, because removing it substantially increases the risk of planet over-subtraction. In the bottom row of Extended Data Fig. 3, we also see that a single OPD plane is not enough model flexibility: in our formulation, the interaction between both OPD planes (separated by a focal plane mask) is important. This illustrates the improvements of FLIP relative to previous autodiff-related efforts, notably ref. ²⁵ (which used a far less flexible model and no model calibration).

Implementation details. We use the Adam optimizer²⁹ with a cosine annealing learning rate scheduler. We use an L1 loss function for our optimization.

We select the learning rates in our model by hand, increasing it across orders of magnitude until we find one that produces a stable optimization with convergence below around 15,000 iterations. Owing to the different magnitudes of different parameters in our model, we assign a different learning rate to each family of parameters. A simple coarse search clearly shows the appropriate scale for each parameter. Extended Data Fig. 4 shows that order-unity variations to learning rates do not substantially affect our results and that increasing or decreasing them by orders of magnitude clearly shows the correct scale for stability and relatively fast convergence. A summary of the order of magnitude of these learning rates per parameter is shown in Extended Data Table 1.

Across different observations, we find that globally scaling all learning rates by some order-unity constant can help with stability and convergence speed. This is probably because of the magnitude of the loss function depending on the brightness of the stars. We find that a global scaling of four works well in all of the experiments presented here, but preliminary analysis on fainter stars showed better results by using 1 or 0.5 instead in those cases.

Similarly, during the last sub-stage of our optimization, we reduce the global scaling factor to around 1% of its original value to further reduce risk of overfitting. Our results are not strongly sensitive to the exact choice of this value but it can be tuned differently for different observations depending on the similarity between reference star and science target.

To determine a stopping point for our optimization, we inject fake planets in the science observation and track their SNR as a function of iteration number. When measuring sensitivity, we inject planets at 5–7 separations and 1–3 planets per separation. For each injection, we inject the fake planet by simulating an off-axis PSF with the nominal STPSF model at the desired location in the field of view and adding it to the science image. We inject them at 4–20 times the 5σ fundamental noise contrast limit. This number may be adjusted depending on the amount of signal loss observed. We inject each planet individually per FLIP run. We measure the SNR every approximately 300 iterations. As expected, the SNR of each injected planet increases for several iterations and then decreases owing to algorithm overfitting. Therefore we stop our optimization at the approximate iteration when the injected planets reach an SNR plateau or peak on average, which was found to be consistent for planets injected at similar separations. Extended Data Fig. 5 shows this for fake planets injected at a variety of separations. In light of this, in our experiments, for simplicity, we choose our stopping point tracking a single fake planet injection. Just as with KLIP, the optimal fake planet injection strategy depends on context. For instance, we can choose as the stopping point the iteration number with maximum average SNR across all planet injections (for example, if the goal is to maximize sensitivity across the field of view) or based on one or several injections at a specific brightness and position (for example, if the goal is to optimize the SNR of a specific planet signal). Extended Data Fig. 5 also shows that our PSF subtraction is not strongly sensitive to our specific choices of fake planet injections.

Experiments

We carried out experiments on four datasets: three with known companions (HIP 65426 in the F444W filter and 51 Eridani in the F410M and F430M filters) and one with no companions, part of a planet search survey (TYC 5899) in the F444W filter. HIP 65426 and TYC 5899 are observed with the MASK335R coronagraph mask, whereas 51 Eridani is observed with the LWBAR coronagraph mask in the NARROW position offset. We selected this diverse set of data to be representative of a wide range of observations with JWST NIRCam coronagraphy. In particular, our four datasets include:

1. 1.

   Three different physical systems (HIP 65426, 51 Eridani and TYC 5899).

2. 2.

   Three different observational filters (F444W, F410M and F430M).

3. 3.

   Two different coronagraphs (NIRCam 335R round coronagraph and NIRCam LWBAR bar coronagraph).

4. 4.

   Both real planets (HIP 65416 b and 51 Eridani b) and simulated planets (TYC 5899).

5. 5.

   Observations with a dedicated reference star with a 9-POINT-CIRCLE dither pattern (HIP 65426), with a dedicated reference star with 5-POINT-BAR dither pattern (51 Eridani) and without dedicated reference stars but with a PSF library composed from a self-referenced survey (TYC 5899). This is relevant for the KLIP performance under a variety of reference star configurations.

6. 6.

   Observations with reference stars close in time to the science target (HIP 65426 and 51 Eridani) and far in time from the science target (TYC 5899).

7. 7.

   Observations of host stars of various spectral types: A2V (HIP 65426), F0IV (51 Eridani) and M3 (TYC 5899).

We run standard PSF subtraction method KLIP and our technique FLIP on all four datasets and compare the performance of each, both with the SNR of the known recovered companions and with the planet detection sensitivity at several separations from the star. We perform our KLIP subtractions in accordance with the optimal set of KLIP hyperparameters described in the respective publications to replicate what the authors considered the best performance. In particular, all KLIP subtractions are done using RDI + ADI (see, for example, ref. ³¹) on both roll angle observations of the science targets, which we then de-rotate and mean-combine. In all cases, we obtain PSF subtractions comparable with the published works. For FLIP, for simplicity, we perform it on only one roll per dataset. Future work could explore the differences in sensitivity between using FLIP only on one roll, FLIP on both rolls and mean-combining the de-rotated images and FLIP on both rolls, using ADI on the resulting images and mean-combining the FLIP + ADI results.

For TYC 5899, the original publication³² presents two epochs: one with a single roll angle (the second failed) and one with two roll angles observed four months after the reference star sequences. We use the latter in this analysis for both KLIP and FLIP, to illustrate the gains in sensitivity when there is substantial drift between reference and science observations. Although ref. ³² used all three observations for their contrast curves and we replicate it only with the last two, we obtain comparable contrast performance: about 3 × 10⁻⁶ at 2.5 arcsec and about 4 × 10⁻⁵ at 1.0 arcsec for KLIP and FLIP and about 3 × 10⁻⁴ for KLIP and about 5 × 10⁻⁴ for FLIP at 0.5 arcsec. Because the target has no known companions, we inject planets at various separations and brightnesses to compare the performance of KLIP and FLIP in a controlled setting.

Extended Data Table 2 shows the optimization hyperparameters used during our reductions with FLIP.

Data preparation

The observations of HIP 65426, 51 Eridani and TYC 5899 and their reference stars were processed with the JWST pipeline (version 1.18.1) and the community-developed tool spaceKLIP (version 2.2.1). Data processing steps (which include ramp-fitting, converting from detector counts to physical flux units and bad pixel cleaning) and observation set-up are detailed in the respective papers for those observations^26,31,32.

The JWST and spaceKLIP pipelines include bad pixel flagging and cleaning. However, imperfect bad pixel cleaning can cause more difficulties in our model-based PSF subtraction than in methods such as KLIP. This is because hot or dead pixels, or variations in pixel responses, cannot be ‘cancelled out’ by subtracting observations with the same pixel variations. Therefore, bad pixels that persist even after the automated cleaning routines are identified by eye and masked during our optimization process.

Choice of reference stars

Both HIP 65426 and 51 Eridani were observations carried out with dedicated reference stars: a bright star with probably no planets around it, of similar spectral type to the science target, observed just before or after the science observation. For these targets, we use as our calibration star (optimization stage 1; Extended Data Fig. 2) the same reference stars used for those observations: HIP 68245 for HIP 65426 and HD 30562 for 51 Eridani. We use a single dither for model calibration with FLIP for simplicity but all dithers are used to build the principal components for KLIP. See refs. ^26,31 for details.

For TYC 5899, the observation was taken as part of a self-referenced programme, in which only science observations are taken, and a large number of these are put together to be used as a reference library for every other science observation. For our KLIP subtractions, we use the full PSF library used in ref. ³² to reproduce the result in their published work. For an apples-to-apples comparison of FLIP and KLIP, we chose one of the observations in the self-referenced library as our calibration star for stage 1 of FLIP. We select the HIP 17695 observation as it has the most similar spectral type to TYC 5899 from the targets in the library. See ref. ³² for details.

PSF subtraction performance metrics

The contrast curves for both KLIP and FLIP were calculated by dividing the image into small concentric annuli and computing the standard deviation within each annulus accounting for small sample statistics⁵⁰ and for planet signal loss owing to both algorithm overfitting and coronagraph throughput. For reference, the signal loss from algorithm overfitting near the location of the measured planets (in the same order as in Fig. 2) for KLIP is about 5–10%, about 20–70%, about 50–80% and about 10–40% and for FLIP it is about 5–20%, about 20–40%, about 30–50% and about 40–60%. We also smooth all datasets with a Gaussian kernel (σ = 1 pixel) before computing the standard deviations, as FLIP residuals are mostly dominated by pixel-to-pixel noise. For HIP 65426, the FLIP contrast curve excludes the central 0.6 arcsec of the image given that a bright lobe of the planet PSF prevents sampling this region, whereas KLIP uses the KLIP-FM algorithm to subtract the planet before computing sensitivity. Future work incorporating planet modelling and subtraction with FLIP will overcome this limitation. See refs. ^26,31,32 for details on computing and calibrating these contrast curves.

For the fundamental sensitivity curve, we use the ERR extensions of the JWST data products, which track the combined total uncertainty in each pixel value from read noise, Poisson noise and flat noise. We generate a random Gaussian sample with mean 0 and σ equal to the total noise level in each pixel. That is, when generating our noise realization, we are simulating what the image would look like after perfect PSF subtraction. We finally compute the contrast curve on this image in the same way as described above for KLIP and FLIP subtractions, assuming no planet signal loss. This assumes that the only source of fundamental noise comes from the science observation. This is essentially the case 0 presented in ref. ⁴. In KLIP, further noise also comes from the reference observations (see ref. ⁴). FLIP is model-based PSF subtraction, so our FLIP subtraction adds no extra fundamental noise. The fundamental noise curves can also be interpreted as a lower limit, or the most optimistic case, of the possible sensitivity. Because we perform FLIP on only one roll per dataset, the noise limit corresponds to the noise limit of that roll.

To supplement Fig. 2, we compute SNR maps for each PSF-subtracted image (Extended Data Fig. 6). We generate a non-coronagraphic PSF model cropped to the central 5 × 5 pixels and cross-correlate our image with the PSF model. Then we compute the SNR of each individual pixel within an annulus of width approximately 0.2 arcsec in a two-step process: we compute the noise per annulus to identify pixels with SNR > 5 as a first pass and then compute a refined noise estimate excluding those pixels. The final SNR map is simply the pixel values of the cross-correlation map divided by the refined noise estimate.

Estimation of the stability-limited sensitivity

In Fig. 3, we show a fiducial stability-limited sensitivity. This limit is computed by performing KLIP on simulated data, in which we can isolate the physical factors that contribute to the KLIP performance limitations. We perform 50 realizations of each observation sequence, stochastically varying the star–coronagraph misalignment and the OPD drift consistent with empirical JWST characterization in flight^5,33, as detailed below. To isolate the effects of these instabilities, we include no other source of noise (for example, photon noise, flat noise) in our simulations.

For the science target, we generate STPSF models using the spectral type of the star and scale the flux within the central 20-pixel radius to the flux of the real data. We simulate two rolls offset by 10° in position angle. For the reference stars, for HIP 65426, we simulate nine PSFs following the 9-POINT-CIRCLE dither strategy, for 51 Eridani five PSFs following the 5-POINT-BAR dither strategy (https://jwst-docs.stsci.edu/jwst-near-infrared-camera/nircam-operations/nircam-dithers-and-mosaics/nircam-subpixel-dithers/nircam-small-grid-dithers#gsc.tab=0) and for TYC 5899, we simulate the same number of reference stars that were used in the PSF library of ref. ³².

We assume that all reference stars are of identical spectral type and brightness to the science target and assume no photon, detector or flat noise. In some targets, this results in the KLIP stability limit being deeper than the fundamental noise limit at large separations from the star (Fig. 3). This is acceptable for our purposes, because at these close separations, we are generally limited by PSF subtraction inaccuracies rather than fundamental noise and we are mostly interested in improving sensitivity in the regions of the image in which we are stability-limited rather than fundamental-noise-limited.

For each observation, we simulate a target acquisition error with a random star–coronagraph misalignment drawn from \({\mathcal{N}}(\mu =0,\) \(\sigma =10)\,{\rm{mas}}\), consistent with target acquisition accuracy in flight³³. We simulate the PSFs with the closest-in-time OPD measurement. We assume that there is no notable OPD drift between reference star dithers. To mimic small OPD drifts owing to telescope slews, for each observation in the real data that required a telescope slew (for example, going from reference star to the science target), we add an OPD drift term to the OPD measurement. The shape of the OPD drift is taken from subtracting the two closest-in-time OPD measurements to the observation being simulated. The root mean square (r.m.s., that is, the amplitude) of the drift term is randomly drawn from \({\mathcal{N}}(\mu =4,\sigma =1)\,{\rm{nm}}\), smaller than drifts measured over day timescales (about 7–9 nm r.m.s.) and bigger than drifts on minutes timescales (about 1–2 nm r.m.s.)⁵.

For HIP 65426, the bright companion has a non-negligible effect on PSF subtraction by skewing the KLIP algorithm (for example, during mean subtraction and in the calculation of Karhunen–Loève (KL) modes). We recreate this effect by injecting a fake planet in our HIP 65426 observation at the same contrast, separation and position angle as the real companion in our simulated observation.

We perform KLIP with ADI + RDI on all targets, using all KL modes for HIP 65426 and 51 Eridani and six KL modes for TYC 5899, following the respective publications. We compute our calibrated stability-limited KLIP sensitivity with injection-recovery tests, injecting planets at 3–7 separations and 2–3 position angles across the 80 × 80-pixel field of view, with the injected flux being four times the uncalibrated KLIP sensitivity at the injection separation. The central 68% range of the resulting contrast curves are plotted as shaded regions in Fig. 3.

TYC 5899 synthetic injections

To compare the sensitivity of FLIP and KLIP on TYC 5899, we inject fake planets with known flux across the image and assess the recovery PSF subtraction. Extended Data Figs. 7 and 8 show a range of planet injections, complementing Fig. 2. At each separation, we inject a planet with flux equal to 20 times the 5σ fundamental noise limit. For PSF subtraction, we process all observations in Extended Data Figs. 7 and 8 identically, using the same learning rates and numbers of iterations for FLIP and the same numbers of annuli and KL modes for KLIP. This is a standard method to quantify sensitivity in KLIP (for example, ref. ³⁷): once a set of annuli and KL modes have been chosen (for example, to recover the signal in Fig. 2), we inject fake planets across the field of view while keeping the annuli and KL modes fixed. This estimates the detection sensitivity across the field of view, given a set of hyperparameters. In FLIP, we use an analogous strategy, by injecting fake planets across the field of view, while keeping the learning rates and numbers of iterations fixed.

We find that FLIP outperforms the standard method KLIP in most cases, especially at closer separations, producing higher SNR on a diverse range of possible planets. We also see that, in this particular dataset, at larger separations, FLIP produces higher signal loss than KLIP in some cases, which explains the contrast curve for KLIP being deeper than FLIP at larger separations in Fig. 3. Future work could address these differences in signal loss. Nevertheless, at these larger separations, we are often background-limited rather than contrast-limited; that is, sensitivity is mainly limited by exposure time rather than by starlight subtraction. In the interesting regime in which we are limited by telescope stability and imperfect starlight subtraction, at smaller angular separations (which also correspond to more Solar-System-like physical separations), FLIP robustly outperforms KLIP.