Aerosols and hydrocarbons in the atmosphere of a white dwarf planet

Data 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 trends^13,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 model³¹. The stellar limb darkening was modelled with the procedures from ref. ³² and a quadratic function using ExoTiC-LD³³. 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. ³⁴. 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 extraction³⁵ 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 model³¹, applying a quadratic limb-darkening law with coefficients generated with ExoTiC-LD³³ 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 methods³⁶ 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 system³⁴. 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.

Grazing transit spectroscopy

The unique transit geometry of the WD 1856 system required us to develop a new approach to express and model transmission spectra. A transmission spectrum encodes the wavelength-dependent effective area of a planet relative to its host star. Exoplanet analyses typically take the spectroscopic planet–star radius ratio from light curve fits, R_p(λ)/R_*, and then express the transmission spectrum as \({R}_{{\rm{p}}}{(\lambda )}^{2}/{R}_{\ast }^{2}\). This quantity is equivalent to the transit depth for a planet with radius R_p(λ) fully occulting a non-limb-darkened star of radius R_*. However, because WD 1856 b is seven times larger than its white dwarf host with a grazing transit, the transmission spectrum cannot be written as \({R}_{{\rm{p}}}{(\lambda )}^{2}/{R}_{\ast }^{2}\). Indeed, the transit depth of WD 1856 b is dependent on time throughout the transit (see ref. ¹¹), with a maximum transit depth corresponding to the time of greatest areal overlap between the planet and its host (Fig. 1). We therefore express the transmission spectrum of WD 1856 b as the wavelength-dependent maximum transit depth at the time of mid-transit, A_p/A_*.

We convert the spectroscopic planet–host radius ratio into the mid-transit transit depth by calculating the time-dependent area overlap of two discs. The overlapping area of two circles with radii R_p (representing the planet) and R_* (representing the white dwarf), separated by a distance d, is given by:

in which

Considering the time of mid-transit (when d = bR_*, in which b is the transit impact parameter), we can express the maximum transit depth as:

in which

We use equations (4), (5) and (6) to map the spectroscopic radius ratio, R_p/R_*, and impact parameter from the spectroscopic light curve fits of each data reduction into the equivalent mid-transit transmission spectrum. We use the uncertainties Python package to propagate errors using these formulae. This approach automatically removes offsets between the different reductions for R_p/R_*, as each corresponding pair of R_p/R_* and b must yield consistent A_p/A_* to have the same transit shape (that is, to match Fig. 1).

We show our final transmission spectra of WD 1856 b, expressed as the mid-transit transmission spectrum (A_p/A_*), in Extended Data Fig. 2. Both reductions clearly detect the strong signature of nightside contamination (the slope to lower transit depths at longer wavelengths) and lead to consistent atmospheric inferences from our retrieval analysis (see the ‘Atmospheric retrieval analysis’ section). We note that the two reductions partially deviate at wavelengths longer than 5 μm—mainly because of differences in the light-darkening treatments and uncertainties in the red edge detector behaviour—so we restricted our atmospheric analysis for WD 1856 b to the NIRSpec PRISM data from 0.5–5 μm.

White dwarf host spectrum

We extracted a calibrated out-of-transit NIRSpec PRISM stellar spectrum for WD 1856 using the FIREFLy data reduction. Starting from the cleaned 2D images, we further flat-fielded, flux-calibrated and extracted the resulting host spectrum. The resulting stellar spectra are shown in Extended Data Fig. 3.

We determined the atmospheric parameters of the host white dwarf by fitting the out-of-transit system flux. We minimized the χ² for a model suitable for cool white dwarfs¹⁴ defined by three parameters: the effective temperature of the white dwarf, T_eff, its photospheric hydrogen-to-helium abundance ratio and the solid angle πR_*²/D². Because the distance D is known from the Gaia DR3 parallax, the solid angle directly constrains the radius of the white dwarf. The radius, in turn, determines the mass and surface gravity of the white dwarf given theoretical white dwarf structure models³⁷. The best-fit solution (Extended Data Fig. 3) corresponds to T_eff = 4,920 K, logg = 8.05 and N_H/N_He = 4.1. This solution yields an Hα line that extends 2% below the continuum, which is consistent with previously obtained optical spectroscopy¹¹ (not considered here in our fit). We also attempted to fit the PRISM spectrum using pure-hydrogen models (that is, without considering N_H/N_He as a free parameter), but the best-fit solution yields a much worse fit to the PRISM data than the mixed H and He atmosphere solution. We calculated limb-darkening coefficients for the best-fitting white dwarf model (using the approach from ref. ³⁸), which we then fixed during the WD 1856 b spectroscopic light curve fits for our two data reductions.

Transmission spectrum modelling

The grazing transit geometry of WD 1856 b, coupled with the clear presence of planetary nightside thermal emission, requires a new modelling approach. The transmission spectrum for a planet with a grazing transit and nightside thermal emission can be written as³⁹:

in which A_p(top) is the area of the planet overlapping the star at the top of the modelled atmosphere (given by equation (1) with R_p = R_p,top), \({{\mathcal{T}}}_{\lambda }\) is the atmospheric transmissivity (\({{\rm{e}}}^{-{\tau }_{\lambda }}\), in which τ_λ is the slant optical depth) in the area element dA and F_p(night),λ and F_*,λ are the observed fluxes from the planetary nightside and white dwarf at Earth, respectively. We reduce the area integral in the first term to a single integral over the fractional annuli of the planet overlapping the star (that is, dA = A_p(r_i,up) − A_p(r_i,low), in which we use the radii of the upper and lower boundaries of each atmospheric layer in place of R_p in equation (1)). The first term in equation (7) represents the wavelength-dependent effective area of the fraction of the planet overlapping the white dwarf, relative to the projected disk area of the white dwarf. The second term accounts for the ‘nightside pollution’/dilution of the transit depth¹⁷ owing to thermal emission from the planetary hemisphere facing the observer.

Transmission spectra of WD 1856 b can also be expressed in terms of emergent fluxes by using the solid angle relation between the observed flux and the emergent (surface) flux, such that equation (7) becomes:

in which R_p,(night),λ is the radius of the emitting thermal photosphere on the nightside (nominally the τ_v,λ = 2/3 pressure level, in which τ_v,λ is the vertical optical depth integrated downwards from the top of the atmosphere). The planet–star surface flux ratio featured in equation (8) is a standard output from radiative transfer codes used to calculate exoplanet emission spectra. Similarly, the transmissivity, \({{\mathcal{T}}}_{\lambda }\), is already calculated by radiative transfer codes calculating standard transmission spectra. Therefore, to calculate the transmission spectra of WD 1856 b, we can construct a model atmosphere and then calculate both the transmissivity from the slant optical depth and the emergent planet–host flux ratio. The observed transmission spectrum then represents a product between a grazing transit transmission spectrum and an ‘upside-down’ emission spectrum.

Atmospheric retrieval analysis

We infer the atmospheric properties of WD 1856 b using the open-source Bayesian atmospheric retrieval code POSEIDON^15,16. We model the atmosphere of WD 1856 b using 100 layers spaced uniformly in log-pressure from 10⁻⁷ to 100 bar. We assume that the atmosphere is well mixed, with consistent atmospheric properties at the day–night terminator and at the nightside, such that only a single set of parameters describe the atmospheric state. We fit for the planetary radius at the 10 bar pressure level and the planetary mass, while fixing the white dwarf radius to R_* = 0.0131 R_Sun and the transit impact parameter to 7.430234 (as the impact parameter uncertainty is already marginalized into the R_p/R_* uncertainties, it does not need to be an independent free parameter). We include the log₁₀ mixing ratios of the following molecules as free parameters: CH₄, NH₃, H₂O, CO₂, CO, HCN, C₂H₂, C₂H₄, C₂H₆, H₂S and PH₃. The remainder of the atmosphere is composed of H₂ and He with an abundance ratio of He/H₂ = 0.17, consistent with the giant planets in the Solar System⁴⁰. We parameterize the temperature profile of WD 1856 b using an adaptation of a prescription used for brown dwarfs⁴¹. This prescription retrieves the temperature at nine pressure nodes (spaced uniformly per decade in pressure from 10⁻⁶ to 100 bar) and interpolates between them with a spline. The nine free parameters defining this temperature profile are the 100 mbar temperature and eight ΔT_i parameters encoding the temperature difference between each pair of nodes. Given the low external irradiation of WD 1856 b, we restrict ΔT_i > 0 to consider physically plausible profiles with temperature monotonically increasing with pressure. Finally, we fit for a three-parameter aerosol model consisting of a power law scattering slope (with exponent γ) and an optically thick cloud-top pressure¹⁵. We do not consider inhomogenous clouds around the terminator, as only a small fraction of the terminator of WD 1856 b occults the surface of the white dwarf during transit.

Our retrieval model is thus defined by 25 free parameters, which we fit using MultiNest’s^42,43,44 Python wrapper PyMultiNest⁴⁵ with 1,000 live points. The priors for each parameter are summarized in Extended Data Table 2. We calculate Bayes factors (that is, odds ratios; \({\mathcal{B}}\)) using Bayesian model comparisons between nested retrieval models, with the retrieval model statistics summarized in Extended Data Table 3. For consistency with the exoplanet literature, we also convert the Bayes factor between two nested models (for example, our reference model and a model excluding CH₄) into an ‘equivalent detection significance’, Nσ, using a standard relation⁴⁶. We note, however, that there are several caveats associated with the Bayes factor to detection significance mapping, so our preferred statistic for model preference is the Bayes factor/odds ratio (see ref. ⁴⁷).

We calculate model transmission spectra of WD 1856 b by solving equation (8) on a wavelength grid ranging from 0.5 to 5.6 μm at R = 20,000. We sample high-resolution pre-computed cross-sections³⁹ onto this wavelength grid, using the following line list sources: CH₄ (ref. ⁴⁸), NH₃ (ref. ⁴⁹), H₂O (ref. ⁵⁰), CO₂ (ref. ⁵¹), CO (ref. ⁵²), HCN (ref. ⁵³), C₂H₂ (ref. ⁵⁴), C₂H₄ (ref. ⁵⁵), C₂H₆ (ref. ⁵⁵), H₂S (ref. ⁵⁶) and PH₃ (ref. ⁵⁷). We also include continuum opacity from H₂ and He collision-induced absorption⁵⁸ and Rayleigh scattering. For the host flux, we use the best-fit white dwarf model shown in Extended Data Fig. 3. Our model transmission spectra are finally convolved with the NIRSpec PRISM point spread function and binned down to the resolution of the observations to calculate the likelihood of each location in the retrieval model parameter space.

Although Fig. 3 compares several retrieved atmospheric properties between the FIREFLy and Juniper data reductions, we provide the full posterior distributions in Extended Data Fig. 4. We find excellent agreement between FIREFLy and Juniper for all retrieved parameters.

To interpret the thermal history of WD 1856 b, we also calculate posterior distributions for the planetary effective temperature from our retrieval results. We calculated the emergent planetary surface flux of WD 1856 b for each atmosphere in the full set of posterior samples from both the FIREFLy and Juniper reductions on a wavelength grid from 1 to 50 μm. For each set of atmospheric parameters, we calculate the corresponding effective temperature using the Stefan–Boltzmann law: \({T}_{{\rm{eff}}}={\left(\frac{1}{{\sigma }_{{\rm{SB}}}}\int {F}_{{\rm{p}},{\rm{surf}},\lambda }{\rm{d}}\lambda \right)}^{1/4}\). Extended Data Fig. 5 shows our retrieved surface flux spectrum for WD 1856 b for both data reductions and the corresponding T_eff posterior distributions. Using the lowest 1σ credible interval (from FIREFLy) and the highest 1σ credible interval (from Juniper), we find a range of 390–412 K for T_eff. We similarly report the 1σ range encompassing both data reductions for M_p in the main text. The roughly 10 K 1σ uncertainty on T_eff is driven by the numerous CH₄ bands detected in our NIRSpec PRISM data setting the relative amplitude of other CH₄ features at longer wavelengths. However, the potential presence of other hydrocarbons, such as C₂H₆, allows a larger surface flux uncertainty in which these species absorb in the mid-infrared (for example, 10–15 μm), which increases the uncertainty in the integrated power and hence T_eff. Longer wavelength observations of WD 1856 b with MIRI LRS/MRS, such as those planned in JWST Cycle 4 (GO-9033 and GO-9157), will constrain T_eff even further.

Mie scattering retrievals

We have established that models including aerosol opacity are required to explain the transmission spectrum of WD 1856 b. Specifically, our free retrieval analysis above infers an opaque cloud deck near 100 mbar and a haze to explain the power-law scattering slope shortwards of 1 μm. The enhanced scattering slope indicates a collection of small particles in the upper atmosphere, but our parametric description is agnostic to the specific aerosol composition. Here we consider retrievals including Mie scattering to investigate which specific aerosol species are consistent with the transmission spectrum of WD 1856 b.

The composition of small, Mie scattering particles can be potentially identified by means of aerosol absorption features at infrared wavelengths, whereas their particle size is encoded by the scattering slope. We assess here which aerosol species and particle sizes can explain the observed scattering slope by means of retrievals including compositionally specific Mie scattering aerosols. We do not test directly for specific species causing the opaque cloud deck, as this deck is probably composed of large particles with muted resonance features^59,60. Because such a condensate cloud deck has no spectroscopic features, it is not possible to determine the composition unless condensates are lofted above the deck and become smaller in size.

We use the Mie scattering retrieval module and database introduced in POSEIDON v1.2 (ref. ⁶¹). Our Mie scattering retrievals use aerosol extinction cross-sections pre-computed from refractive indices. We mainly consider a simple aerosol model parameterized by the log₁₀ mean particle size (logr_m \({\mathcal{U}}\,[-3,1]\)) and aerosol log₁₀ volume mixing ratio (log aerosol \({\mathcal{U}}\,[-30,-1]\))—representing a well-mixed aerosol uniformly distributed within the atmosphere. We also tested more complex aerosol models that fit for pressure-dependent aerosol mixing ratios but these all reduced to a pressure-independent model. Our Mie retrievals use a six-parameter pressure–temperature (P–T) profile⁶². We conduct these Mie retrievals on the FIREFLy data reduction.

We ran retrievals with a suite of aerosol species representing three different aerosol formation regimes that could be relevant in upper atmosphere of WD 1856 b. The first aerosol regime represents disequilibrium hazes and soot species that can be produced by photochemistry: Titan tholins (tholins^63,64), carbon soot (C (ref. ⁶⁵)), water-rich organic haze at two temperatures (ExoHaze 300K, ExoHaze 400K (ref. ⁶⁶)) and hexene (C₆H₁₂ (ref. ⁶⁷)). The second aerosol regime represents the myriad of sulfide and chloride clouds that form in brown dwarfs at the T–Y transition (400–1,300 K)²⁷ alongside Cr: chromium (Cr; Lynch and Hunter in ref. ⁶⁸), magnesium sulfide (MnS (ref. ⁶⁹)), sodium sulfide (Na₂S (ref. ²⁷)), zinc sulfide (ZnS (ref. ⁷⁰)) and potassium chloride (KCL; Palik and Addamiano in ref. ⁷¹) (ordered by condensation temperature). The third aerosol regime consists of condensed ices that form deep cloud decks in Solar System planets and potentially cooler Y dwarfs (≤400 K)^72,73. These ices could cause the opaque cloud deck found in our retrievals above, which are then lofted to higher atmospheric pressures to cause the observed scattering slope, or they could condense in situ in the colder upper atmosphere: water ice (H₂O (ref. ⁷⁴)), ammonia ice (NH₃ (ref. ⁷⁵)) and methane ice (CH₄ (ref. ⁷⁶)) (ordered by condensation temperature).

We find that all aerosol species, with the exception of MnS and hexene, provide good fits to the scattering slope and only imprint weak absorption features into the transmission spectrum. Using Bayesian model comparisons, the best-fit haze and soot species is the water-rich organic ExoHaze (the 400K variant), the best-fit T–Y dwarf cloud species is KCl and the best-fit ice is NH₃. Of these three aerosols, KCl has the highest Bayesian evidence. The potential presence of KCl would be consistent with expectations for cold T–Y dwarf models²⁷, in which KCl forms the highest, low-density cloud. However, we note that a simple grey cloud deck + haze model (as used in the main text) is preferred over KCl by about 2σ. Therefore, the present data for WD 1856 b is not sufficiently precise to identify a clear preference for which specific aerosols are present in the atmosphere of WD 1856 b.

Our Mie scattering retrievals provide insights into the range of particle sizes and abundances compatible with the short wavelength scattering slope WD 1856 b (Extended Data Fig. 6). The ExoHaze and NH₃ ice models favour a collection of small particles (about 0.03 μm) with low mixing ratios (about 10⁻¹⁴), whereas the KCl model favours even smaller particles (roughly 0.01 μm) with a higher abundance (about 10⁻⁸). Compared with our default grey cloud deck + haze retrieval model, we find consistent results for other model parameters to within 1σ. In particular, we show that the retrieved planetary mass is not sensitive to the assumed aerosol model. We do find roughly 1 dex lower median CH₄ abundances for the Mie scattering retrievals and hence a lower C/H ratio, but the CH₄ abundance distribution is still consistent with our results in the main text. We note that the marginal evidence of C₂H₆ strengthens when including Mie scattering compared with the deck + haze model (Extended Data Fig. 6) but this molecule is not strongly detected with the present data.

Our retrieved temperature structure from the Mie scattering retrievals also indicates an atmosphere that is much warmer than the equilibrium temperature of WD 1856 b (Extended Data Fig. 6). As with our grey cloud and haze retrieval, we also find a temperature of about 400 K in the thermal photosphere near 10–100 mbar. However, because the Mie scattering retrievals cannot produce an optically thick cloud deck at the pressures required to obscure thermal emission from the deep atmosphere (approximately 10^−1.5 bar), the Mie retrievals compensate by making the P–T profile essentially isothermal in the deep atmosphere (that is, the retrieved P–T profile shown in Fig. 2 is more physical). We note that our uniform aerosol Mie scattering retrievals are incompatible with the P–T profile used in the main text⁴¹, as a collection of small aerosols are not able to simultaneously block the deep adiabatic thermal flux and fit the scattering slope. The P–T profile parameterization chosen here⁶² tends to favour a nearly isothermal upper atmosphere, which suffices for the exploration of the aerosol properties consistent with the scattering slope of WD 1856 b. Future explorations of the cloud structure and radiative properties of WD 1856 b, such as composite cloud models with multiple scattering, are a rich area to deepen our understanding of the atmosphere of WD 1856 b.

Evolution of the WD 1856 system

Host progenitor and white dwarf

We examined the evolution of the progenitor star of WD 1856 by consulting the MIST evolutionary models⁷⁷ for non-rotating solar-metallicity stars in the appropriate mass range (\({M}_{{\rm{progenitor}}}=1.3{6}_{-0.18}^{+0.29}\,{M}_{{\rm{Sun}}}\)). From these models, we extracted fiducial estimates of the main sequence lifetime (\({4}_{-1.8}^{+2.4}\,{\rm{Gyr}}\)) using an initial–final mass relation^78,79, the duration of the thermally pulsing AGB stage (\(1.5{5}_{-0.10}^{+0.26}\,{\rm{Myr}}\)) and the post-AGB/pre-white-dwarf stage (\(0.03{4}_{-0.002}^{+0.053}\,{\rm{Myr}}\)). The latter is defined here as the elapsed time between the final thermal pulse and the cooling of the exposed core to an effective temperature of 100,000 K. This yields a total system age of \(9.{4}_{-1.9}^{+2.5}\,{\rm{Gyr}}\).

We calculated the cooling age of the white dwarf host by evolving MESA white dwarf models of the appropriate mass down to T_eff = 4,920 K. We used MESA r23.05.1 (ref. ⁸⁰). This MESA release now includes carbon–oxygen fractionation⁸¹, which is important here as the white dwarf is in the process of crystallizing. A standard helium layer of \(\log {M}_{{\rm{He}}}/{M}_{* }=-2\) was assumed, whereas a relatively thin hydrogen layer of \(\log {M}_{{\rm{H}}}/{M}_{* }=-6\) was used. This is much thinner than the canonical value of \(\log {M}_{{\rm{H}}}/{M}_{* }=-4\) (ref. ⁸²) but is motivated by the fact that the model atmosphere analysis points to an atmosphere containing a mix of hydrogen and helium. This presumably requires the superficial convection zone to extend just below the hydrogen layer, thereby diluting hydrogen with helium. From this constraint, we can estimate \(\log {M}_{{\rm{H}}}/{M}_{* }=-6\) (ref. ⁸³). Cooling calculations were performed for different carbon–oxygen core composition profiles to account for current model uncertainties: a standard profile predicted by stellar evolution⁸⁰ and an asteroseismologically derived stratification⁸⁴ were used. We also calculated cooling models using different electron thermal conductivities^85,86 to account for current uncertainties at the transition between the regimes of moderate and strong degeneracy⁸⁷. From this analysis, we find a cooling age of 5.4 ± 0.7 Gyr, in which the uncertainty includes the systematic uncertainty sources listed above and a 2% uncertainty on the T_eff of the star and mass typical of white dwarfs in this temperature range⁸⁸. This cooling age is consistent with estimates produced by other stellar evolution codes^82,89.

Thermal history of WD 1856 b

We reconstructed the thermal evolution of WD 1856 b under the assumption that the cooling of the planet after migration has been similar to the cooling undergone by a substellar object after formation. We used cooling models from the ATMO2020 (ref. ⁹⁰) and Sonora Bobcat⁹¹ model grids. Each grid tabulates global quantities such as luminosity, effective temperature, radius and surface gravity as a function of age for substellar objects of a given mass and bulk chemical composition, starting from an initial condition with high entropy. Both provide self-consistent evolutionary–atmospheric modelling frameworks, in which the structure and evolution of the fully convective, adiabatic interior are computed with a cloudless, non-grey, rainout-chemical-equilibrium atmosphere as the surface boundary condition. The most important difference between ATMO2020 and Sonora Bobcat is that the former neglects some relevant opacity sources at effective temperatures above 2,000 K, leading to faster cooling at high temperatures in the ATMO2020 models. ATMO2020 provides models of solar-metallicity objects, whereas Sonora Bobcat provides models for both solar-metallicity and metal-enriched ([M/H] = +0.5) objects. We consider these three sets of models in our analysis below.

Reconstructing the thermal history of WD 1856 b requires us to choose a model grid and specify three parameters: planetary mass (M_p), current planetary effective temperature (T_eff,p) and current white dwarf cooling age (t_wd). For a given M_p and model grid, we obtained the effective temperature as a function of time by adding a uniform offset (t₀) to the model age (t_p) such that the model temperature matches T_eff,p at t_p = t_wd − t₀ (using linear interpolation between the tabulated model ages and temperatures). The cooling models predict a high effective temperature (about 1,500–3,000 K) at t₀; these values are plausible for planets that have been tidally heated during high-eccentricity migration^92,93 or have survived a common-envelope phase^8,9,94. We therefore interpret t₀ as an estimate of the time of the planet’s reheating during migration, expressed as a white dwarf cooling age. Because cooling is rapid at high T_eff,p, our estimate of t₀ is robust to theoretical uncertainties in what temperature the planet should be immediately after migration.

We considered cooling models with M_p between 0.5 and 20 M_J, covering the range of samples from the atmospheric retrieval posterior distributions for the FIREFLy and Juniper JWST data reductions. Each grid samples a finite number of mass values; when considering objects of arbitrary mass between grid points, we used the cooling model with the nearest mass on the grid. Both the ATMO2020 and Sonora Bobcat models are spaced by about 0.5–1.0 M_J in mass over the range we consider, so our approach does not introduce substantial error in a given reconstruction compared with interpolating between adjacent models.

We generated ensembles of possible thermal histories using the mass and effective temperature constraints derived from the NIRSpec PRISM transmission spectrum of WD 1856 b. Specifically, we considered nearly 10,000 values of M_p from the atmospheric retrieval posterior distribution alongside the nearly 10,000 corresponding values of T_eff,p obtained using the procedure described above for each data reduction (see the ‘Atmospheric retrieval analysis’ section). Our samples of M_p and T_eff,p are not statistically independent, as each pair of values is derived from a single sample from the distribution of atmospheric models consistent with our NIRSpec PRISM transmission spectrum. This has an important influence on the range of thermal histories that we can infer from the data, because the cooling rate is a sensitive function of mass.

On the other hand, our estimated t_wd is independent of our atmospheric retrieval analysis. For each pair of M_p and T_eff,p values, we generated ten random values of t_wd drawn from a Gaussian distribution with mean 5.4 Gyr and standard deviation 0.7 Gyr. Each ensemble therefore comprises about 100,000 possible thermal histories consistent with the transmission spectrum of WD 1856 b. We generated one ensemble for each model grid. Extended Data Fig. 7 shows the distribution of calculated t₀ values for the three cooling models and two data reductions, from which we derive a statistical constraint on t₀. The results reported in the main text were obtained using the solar-metallicity Sonora Bobcat models (solid orange and green histograms; see also Fig. 4). If we use the ATMO2020 models, we find a comparable \({t}_{0}=4.{3}_{-1.1}^{+0.9}\,{\rm{Gyr}}\) for FIREFLy and \(4.{6}_{-1.0}^{+0.8}\,{\rm{Gyr}}\) for Juniper. Using the metal-enriched Sonora Bobcat models yields \({t}_{0}=4.{2}_{-1.4}^{+1.0}\,{\rm{Gyr}}\) for FIREFLy and \(4.{5}_{-1.1}^{+0.9}\,{\rm{Gyr}}\) for Juniper. The conclusions we draw from modelling the thermal evolution of WD 1856 b are therefore robust to both the JWST data reduction and the choice of cooling models, given the models available at present.

In a small fraction of cases (<0.15%) for each ensemble, we calculate values t₀ < 0. These correspond to the highest M_p values sampled from the atmospheric retrieval posterior. Negative values of t₀ arise in these cases because we have calculated t₀ by extrapolating the cooling models back to the effective temperatures expected among newborn brown dwarfs of about 20 M_J (>2,000 K). However, these results are unphysical according to the interpretation of t₀ as the time elapsed between the end of the AGB phase and the reheating/migration of the planet. If we stipulate that reconstructed thermal histories be truncated for t₀ < 0, then these few cases are consistent with common-envelope evolution in that it is possible for the planet to have achieved its current temperature by passively cooling since the end of the AGB phase (albeit from a cooler, lower-entropy state than those implied in cases with t₀ > 0). The fact remains, however, that most of the cases (>99.85%) imply t₀ values that cannot coincide with a common-envelope phase in all three ensembles. Thus, we conclude that reheating during the white dwarf phase (consistent with high-eccentricity migration) is preferred over reheating during common-envelope evolution at >2σ (for FIREFLy) and >3σ (for Juniper). Further theoretical study is needed to corroborate or qualify this conclusion, as we describe below.

Our method of reconstructing thermal histories is based on backward extrapolation of the effective temperature only. However, the cooling models also predict the evolution of the radius of WD 1856 b; these predictions should agree in principle. In Extended Data Fig. 8, we show the radius evolution implied by our reconstruction method for 100 samples from the Sonora Bobcat ensemble for both data reductions. For the observed radius, we use the best-fitting value R_p = 0.911 ± 0.020 R_J from the FIREFLy reduction. We also include a systematic error of ±0.050 R_J given by the range of best-fitting radius values covered by the two data reductions, for a total uncertainty of ±0.054 R_J. We see that many of the temperature-based reconstructions overestimate the radius of WD 1856 b by about 2σ. Future efforts to understand the thermal evolution of WD 1856 b should reproduce both the effective temperature and radius. A clue as to the origin of this discrepancy comes from the heavy-element enrichment of the envelope of WD 1856 b, suggested by our retrieved CH₄ abundance, as planetary radius decreases with increasing metallicity at a fixed mass and internal entropy. Model grids of comparable quality with ATMO2020 and Sonora Bobcat that are applicable to objects as massive (about 7 M_J) and metal-rich (about 100 × solar) as WD 1856 b have not been developed or published to our knowledge.

We note that we neglected the irradiation of WD 1856 b by the host white dwarf in our reconstruction of the planet’s thermal history. Irradiation is a key ingredient in modelling the structure and evolution of short-period exoplanets around main-sequence stars, such as hot Jupiters^95,96. The importance of irradiation in the case of WD 1856 b can be gauged by calculating the ratio of the power emitted from the photosphere of the planet to the power incident on the planet from the star:

in which T_eff,* and T_eff,p are measured effective temperatures of the host star and planet, respectively, R_* is the host radius and a_orb is the orbital semimajor axis (assuming a near-circular orbit). Using the system parameters as determined in this work, we calculate \({\mathcal{R}}\approx 25\), indicating that the self-luminosity of the planet overwhelms the power received from the star. Our reconstructed histories generally find that this ratio was larger in the past (except perhaps in the first several Myr after the white dwarf formed). Thus, we argue that irradiation has had a small effect on the previous thermal evolution of WD 1856 b. It would be of interest to self-consistently model the evolution of a substellar body with time-dependent irradiation, as would be the case in proximity to a cooling white dwarf. We leave this for future work.

Alternatives to reheating during migration

We considered several alternative explanations for the increased effective temperature of WD 1856 b, all of which we deemed implausible or unlikely. We briefly describe each of them here, along with our reasoning.

First, the observed effective temperature of WD 1856 b cannot be explained purely by passive cooling over the system’s total age of about 10 Gyr (ref. ¹). This is readily ruled out by consulting theoretical cooling models^90,91. To have an effective temperature of about 400 K at an age of 10 Gyr, WD 1856 b would need to have a mass of approximately 24 M_J. Our observations rule out such a high mass at >3σ confidence.

The mass of WD 1856 b may be above the threshold for deuterium fusion in its core (about 13 M_J) within 2σ. However, although it is possible that WD 1856 b was once heated internally by nuclear reactions, this cannot explain its present-day properties. Models of deuterium-burning brown dwarfs predict a total luminosity many orders of magnitude greater than that of WD 1856 b (ref. ⁹⁷). The duration of deuterium burning (about 3–50 Myr depending on mass⁹⁷) is much shorter than the total age of the system, so the primordial deuterium WD 1856 b would have been destroyed early in the main-sequence lifetime of the host.

Owing to the proximity of WD 1856 b to its host, tidal interactions are another possible heat source inside the planet; this would be analogous to the heating of the Galilean satellite Io through its tidal interaction with Jupiter⁹⁸. For tidal heating to operate, the orbit of WD 1856 b would need to be slightly eccentric rather than circular as is typically assumed. Assuming that the power dissipated by tidal friction is equal to the total power emitted by WD 1856 b, we calculate the effective temperature of the planet, using the standard ‘equilibrium tide’ theory⁹⁹, as:

Here G is the gravitational constant, σ_SB is the Stefan–Boltzmann constant, M_* is the mass of the host white dwarf, R_p is the radius of WD 1856 b and a_orb and e_orb are respectively the orbital semimajor axis and eccentricity (with e_orb ≪ 1). The quantities k_2p and τ_p, known respectively as the tidal Love number and tidal lag time, characterize the dissipation inside WD 1856 b in the standard equilibrium tidal theory. The reference values of M_*, R_p and a_orb used in equation (10) match the observed system parameters. For k_2p and τ_p, we use values similar to those inferred for Jupiter’s dissipation of the tide raised by Io^100,101. We see that tidal heating could, in principle, sustain the observed effective temperature of WD 1856 b for an orbital eccentricity of about 0.02 (the highly uncertain values of k_2p and τ_p notwithstanding). This would be consistent with orbital circularization in the end stage of high-eccentricity migration. However, the same dissipation would damp the orbital eccentricity on a characteristic timescale of roughly 0.075 Gyr (ref. ⁹⁹). In this picture, we are observing WD 1856 b in the very last, short-lived stage of high-eccentricity migration. Although we cannot rule it out based on the available data, we consider this explanation unlikely.