GoKawiilData reduction
Our PRISM observation used 33 groups across 240 integrations, with each integration lasting 29.8 s, for a total observing time of 1.98 h. Although the transit duration of WD 1856 b only lasts 8 min, we selected this observing window to ensure JWST captured a transit of WD 1856 b with sufficient out-of-transit baseline for detector settling. We used two independent codes to reduce the NIRSpec PRISM observation of WD 1856 b and extract transmission spectra: FIREFLy and Juniper. Here we detail the reduction and light curve fitting approach used by each code.
FIREFLy
We first reduced the WD 1856 b data using the Fast InfraRed Exoplanet Fitting Lyghtcurve (FIREFLy)13,29,30 reduction suite. The reduction started with the uncalibrated (uncal.fits) images and a customized jwst pipeline reduction. During stages 1 and 2, the 1/f noise was removed at the group level, using the top and bottom six rows to measure the count level for each column and subtracting the median value. The dark current step was skipped and the jump step was performed with a rejection threshold of 20. During stage 3, we then used the custom-run pipeline 2D images after the jwst.assign_wcs step and performed customized cleaning of bad pixels, cosmic rays and hot pixels. We used cross-correlation to measure the positional shift of the spectral trace across the detector and shift-stabilized the images with flux-conserving interpolation. This procedure has been found to reduce the amplitude of position-dependent systematic trends13,29. An aperture size of 5.7 pixels was used to extract the spectra, with this product used to fit the transit light curves and extract the exoplanet spectra.
During stages 4 and 5, we fit the white light curve using a linear baseline and a limb-darkened transit model31. The stellar limb darkening was modelled with the procedures from ref. 32 and a quadratic function using ExoTiC-LD33. The values were fixed to the best-fitting theoretical host white dwarf model values (see the ‘White dwarf host spectrum’ section below). We measured the white light curve using the 0.55 to 3 μm region, such that the planetary emission would not bias the transit depth and resulting system parameters. The NIRSpec PRISM spectral time series for FIREFLy is shown in Extended Data Fig. 1). Several pre-transit exposures showed abnormally low flux levels, which we flagged as outliers and removed from the remaining analysis. These outliers seem to be because of clusters of bright/hot pixels, so are probably associated with snowball events. We fixed the period using the results from ref. 34. The best-fit white light curve system values are given in Extended Data Table 1. The transmission spectral light curves used 3-pixel binning and were fit with the same model, setting the system parameters to be fixed to the white light curve values.
Juniper
We also applied Juniper, a new custom pipeline for JWST NIRSpec observations, to reduce our WD 1856 b data. Juniper contains a wrapper for stages 1 and 2 of the jwst pipeline with custom steps. We start our processing from Juniper stage 1 with the uncal.fits files from the Mikulski Archive for Space Telescopes (MAST). We opt to disable the jwst stage 1 jump detection step and instead handle cosmic rays through custom procedures in later stages. Before ramp fitting, group-level background subtraction is performed using the top and bottom six rows as the background region to reduce scatter in the extracted light curves. We spatially filter 3σ outliers from this region and average along columns to determine the background level of counts per column, which is subtracted from the full group. We then proceed with jwst stage 1 ramp fitting and gain scaling. Juniper stage 2 is a pure wrapper for jwst stage 2; we carry out this stage with the flat field and photom steps disabled, as neither is required to measure transit depth and we observed the former to increase the noise in the spectral extraction.
Juniper stage 3 performs extra cleaning at the integration level. We first mask pixels flagged by the jwst pipeline for data quality issues. We then reject cosmic rays in time over two iterations, replacing 6.5σ outliers with the median value of the pixel in time. Finally, a second round of background subtraction, using the same strategy as the group-level background subtraction in stage 1, is performed using the top and bottom three rows with outliers masked at 3σ.
Stage 4 of the Juniper pipeline extracts 1D spectra, which are subsequently binned to produce light curves. Our aperture is centred on the brightest row of the trace and extends ±3 pixels above and below it. We perform optimal extraction35 to extract the 1D spectrum, taking as our extraction profile the median image of the trace contained in the aperture, normalized along columns. We sum across all pixels from 0.552 to 3 μm to extract a broadband light curve; we choose not to include light from wavelengths longer than 3 μm as this light is strongly affected by contamination from nightside thermal emission, which affects the determined system parameters (for example, semimajor axis, transit epoch, impact parameter). We then bin every 3 pixels to produce 137 median-normalized spectroscopic light curves at nearly native resolution, spanning 0.552 to 5.360 μm. Our extracted broadband and spectroscopic light curves would typically be sigma-clipped to further remove outliers; however, this technique is prone to clipping out the transit itself owing to the short transit duration and large transit depth. We therefore disable this procedure and use alternative outlier rejection procedures in stage 5.
Juniper stage 5 is the final stage of the pipeline, which fits transit models to each light curve to extract transit depth and produce a transmission spectrum. Our fitting procedure is a two-step process combining linear and nonlinear fitting techniques, which we use to clean outliers that sigma-clipping cannot safely remove. We start by using a linear least-squares estimator to fit a batman transit model31, applying a quadratic limb-darkening law with coefficients generated with ExoTiC-LD33 using a custom white dwarf model produced by fitting a white dwarf spectrum to the out-of-transit flux of WD 1856 (described below). Our model is further multiplied by a linear-in-time trend to account for visit-long ramp effects. We then compute the residuals of the fitted transit and systematics model and clip any points in the light curve that produce 3σ outliers in the residuals. We compute the standard deviation of the sigma-clipped residuals to estimate the photometric uncertainty, which we supply to the next step of our fit procedure. We refit the sigma-clipped light curve using Markov chain Monte Carlo methods36 to extract our final planet–star radius ratio spectrum. We first fit our 0.552 to 3 μm broadband light curve using this two-step fitting process to determine the broadband depth, semimajor axis a/R * , inclination i, mid-transit epoch t C and linear-in-time systematics model parameters, from which we derive the impact parameter b and its uncertainty. The orbit period P was held fixed to 1.407939217 days based on a follow-up paper studying transit timing variations in the WD 1856+534 system34. Our broadband light curve analysis yielded a/R * = 339.25 ± 5.92 and b = 7.34 ± 0.20. We then fix these values as well as t C as determined by the broadband light curve fit for all subsequent spectroscopic light curve fits. We fit every spectroscopic light curve with our two-step process to determine R p (λ)/R * and the linear systematics trend in every wavelength channel. Our broadband light curve fit achieves residuals of 522 ppm, whereas our spectroscopic fits achieve median residuals of 8,153 ppm. We present our fitted system parameters (a/R * , b, R p /R * ) and broadband transit depth in Extended Data Table 1.
... continue reading