Aerosols and hydrocarbons in the atmosphere of a white dwarf planet

Abstract

Most stars, including our Sun, will one day evolve into red giants and, subsequently, white dwarfs. Several planet candidates have recently been identified orbiting white dwarfs^1,2,3,4, demonstrating that planets can survive the stellar post-main-sequence stage intact. Little is known about the atmospheric composition of post-main-sequence planets, with the most evolved transiting planets with atmospheric detections so far orbiting subgiants^5,6. Here we report an atmospheric detection for the white dwarf planet WD 1856 b, achieved through transmission spectroscopy with the James Webb Space Telescope (JWST) Near-Infrared Spectrograph (NIRSpec) PRISM. Our 0.5–5.0-μm spectrum reveals the presence of hydrocarbons (odds ratio of 167:1–5,377:1, with CH₄ preferred at 17:1–30:1), aerosols (2 × 10⁵:1–2 × 10⁶:1) and thermal emission from the planetary nightside (2 × 10⁶³:1–2 × 10⁷³:1). Our spectral analysis constrains the mass of WD 1856 b to 4.3–10.9 M_J, finds a carbon-enriched atmosphere (with a CH₄ abundance of approximately 7%) and an effective temperature exceeding the expected planetary equilibrium temperature (390–412 K versus 160 K). On the basis of cooling models, these results indicate that WD 1856 b underwent a migration-related reheating event 3.0–5.5 Gyr into the white dwarf phase, consistent with post-main-sequence tidal evolution to the present-day 0.02-au circular orbit. Our results provide a window into the ultimate fate of giant planets orbiting stars with masses similar to our Sun.

Main

We observed a transit of WD 1856 b on 27 April 2023 with JWST’s NIRSpec instrument, using the PRISM mode, as part of Guest Observer Program 2358 (PI: Ryan MacDonald). WD 1856 b is a cool (T_eq = 160 K), Jupiter-sized (0.9 R_J) planet transiting the white dwarf WD 1856+534 (ref. ¹). The host white dwarf (T_eff = 4,710 K, 0.0131 R_Sun, 0.518 M_Sun, t_cool ≈ 6 Gyr (ref. ¹)) is 25 pc from Earth and has a DA spectral class⁷. The current close-in orbit of WD 1856 b (0.02 au/1.4 days) requires orbital evolution after the main sequence to avoid engulfment during the red giant phase. Hypotheses for the orbital evolution of WD 1856 b include common-envelope evolution during the red giant or asymptotic giant branch (AGB) phase^8,9 or, alternatively, high-eccentricity migration¹⁰. Distinguishing these scenarios has proved challenging, given the uncertain planetary mass^1,11,12. As the host white dwarf WD 1856 is a relatively faint target (J = 15.677), we gave priority to the wide wavelength range (0.6–5.5 μm at a mean resolving power of λ/Δλ ≈ 100) and high throughput of the NIRSpec PRISM over other JWST instrument modes. Our observation lasted 1.98 h, of which the transit of WD 1856 b lasted 8 min.

We reduced the JWST observations using two independent data pipelines, FIREFLy¹³ and Juniper (newly presented here; Methods). Each reduction pipeline yielded sets of spectroscopic transit light curves for WD 1856 b, which we independently fit with a transit model (Methods). Because most existing stellar model grids do not include models for evolved or remnant objects, our transit models used a custom limb-darkening prescription for WD 1856 derived from a white dwarf model¹⁴ fitted to the out-of-transit host flux (Methods). Figure 1 shows transit light curves of WD 1856 b from FIREFLy, integrated over two broadband wavelength regions from 0.55–1.71 μm and 4.19–4.96 μm. Our JWST observations find a 3% shallower transit depth in the second, near-infrared, wavelength region compared with visible wavelengths. The physical mechanism for this transit depth difference is planetary thermal emission, as we establish below. Our PRISM transmission spectrum of WD 1856 b from the FIREFLy code is shown in Fig. 2a,b.

Fig. 1: Detection of nightside thermal emission from the atmosphere of WD 1856 b.

a, JWST NIRSpec PRISM transit observations of the white dwarf planet WD 1856 b integrated over two broadband wavelength regions: 0.55–1.71 μm (blue data points) and 4.19–4.96 μm (red data points). The best-fitting model transit light curves are overplotted (blue and red curves). The transit event is shallower at longer wavelengths owing to planetary thermal emission diluting the transit depth. The top schematic depicts the grazing transit geometry of the WD 1856 system, with the planet and white dwarf shown to scale. b, Zoom-in of the mid-transit to demonstrate the high signal-to-noise detection of the nightside thermal emission. The 1σ errors for the blue data points are too small to be seen (the mean 1σ errors are 5.3 × 10⁻⁴ from 0.55–1.71 μm and 2.8 × 10⁻³ from 4.19–4.96 μm). c, Residuals between the data and best-fitting model.

Source data

Fig. 2: Atmospheric retrieval of the transmission spectrum of WD 1856 b.

a, The JWST NIRSpec PRISM transmission spectrum of WD 1856 b (FIREFLy data reduction; orange circles with 1σ error bars) is compared with the retrieved model spectrum (purple line and contours, showing the median, 1σ and 2σ credible intervals). Several absorption features from CH₄ are detected, alongside nightside thermal emission, continuum aerosol opacity and tentative evidence of PH₃ and C₂H₆. b, Zoom-in to highlight the short-wavelength spectral features. c, Retrieved temperature profile (purple contours) and cloud-top pressure (grey contours) from the transmission spectrum of WD 1856 b. An optically thick cloud deck near 100 mbar (grey hatching) blocks thermal emission from deeper layers. The top of the hatched region corresponds to the median retrieved cloud-top pressure and the 1σ and 2σ credible regions appear above (horizontal grey shading). Most thermal emission arises from the cloud-top pressure near 100 mbar (red arrows).

Source data

We model the transmission spectrum of WD 1856 b using the radiative transfer and retrieval code POSEIDON^15,16, adapted here to include partial planet–star overlap (that is, grazing transit geometry) and the effect of nightside thermal emission^17,18 (Methods). We consider the main carbon, oxygen, nitrogen, sulfur and phosphorous carriers expected in a H_2-dominated atmosphere at the equilibrium temperature of WD 1856 b and several disequilibrium tracers, resulting in the following gases: CH₄, NH₃, H₂O, CO₂, CO, C₂H₂, H₂S, PH₃, HCN, C₂H₄ and C₂H₆. Our retrievals freely fit the temperature structure, chemical composition, an opaque cloud deck and a power law haze¹⁵. Given the uncertain mass of WD 1856 b, our retrievals fit the mass jointly with the atmospheric properties. We conducted retrievals on both data reductions, finding excellent consistency between FIREFLy and Juniper (Methods). The retrieved spectrum and temperature structure for the FIREFLy reduction are shown in Fig. 2.

Our retrievals establish that at least one hydrocarbon ($4.5\sigma /\mathrm{ln}{\mathcal{B}}\,=$ $8.59$ for Juniper, $3.6\sigma /\mathrm{ln}{\mathcal{B}}=5.12$ for FIREFLy) is present in the atmosphere of WD 1856 b. The main contributor to the hydrocarbon inference is CH₄ ($3.1\sigma /\mathrm{ln}{\mathcal{B}}=3.40$ for Juniper, $2.9\sigma /\mathrm{ln}{\mathcal{B}}=2.87$ for FIREFLy) but C₂H₆ also contributes to the best-fitting model. Figure 3a,c shows spectral contributions to our best-fitting model, which achieves an excellent fit to the FIREFLy data (${\chi }_{\nu }^{2}=1.03$ with 102 degrees of freedom). The fit quality is slightly less good for Juniper (${\chi }_{\nu }^{2}=1.28$, indicating consistency with the data within 2σ for 102 degrees of freedom) but the retrieval results are nonetheless consistent between the two data reductions. CH₄ causes three prominent absorption bands near 1.75 μm, 2.3 μm and 3.3 μm, with the retrieved CH₄ abundance (2–20%) indicating an atmospheric metal enrichment of about 100 × solar ($\log ({\rm{C}}/{\rm{H}})=2.2{4}_{-0.40}^{+0.33}$ for FIREFLy and $2.2{1}_{-0.44}^{+0.35}$ for Juniper in solar units, with corresponding 3σ lower limits of 10.8 × solar and 7.9 × solar, respectively). We confirmed that the CH₄ inference arises from several bands by repeating a FIREFLy retrieval with the data from 3.1–3.5 μm masked, yielding a 3.0σ model preference. The best-fitting model includes a contribution from PH₃ near 4.3 μm, which could be an indicator of disequilibrium vertical mixing as in the upper atmosphere of Jupiter¹⁹. However, model degeneracies and lower signal-to-noise at longer wavelengths preclude a detection of PH₃ with the present data. Similarly, a feature attributable to C₂H₆ at 3.4 μm is partially degenerate with CH₄ absorption near 3.3 μm, so C₂H₆ is not significantly detected. However, the combined evidence for any hydrocarbon (that is, CH₄, C₂H₂, C₂H₄ and C₂H₆) is statistically significant at approximately 4σ. No other chemical species (for example, NH₃ or H₂O) are detected (see Methods and Extended Data Fig. 4 for upper limits).

Fig. 3: Best-fitting spectrum and atmospheric properties of WD 1856 b.

a, Spectral decomposition illustrating the model components required to explain the JWST transmission spectrum of WD 1856 b (FIREFLy reduction; orange circles with 1σ error bars). The best-fitting model (black curve) contains absorption from several CH₄ features (purple shading), H₂ + He continuum absorption (dark grey shading) and tentative evidence of C₂H₆ near 3.4 μm (brown shading) and PH₃ near 4.3 μm (green shading). Aerosol continuum opacity (light grey shading) blocks thermal emission from the deep atmosphere in windows between CH₄ bands, obscuring the ‘negative spikes’ that would otherwise be observed near 0.7 μm, 2.9  μm and longwards of 4.0 μm for a clear atmosphere. The best-fitting model without nightside thermal emission (dashed blue curve) cannot explain the much lower transit depths longwards of 3.5 μm. b, Retrieved atmospheric properties, planetary mass and aerosol properties from the transmission spectrum of WD 1856 b for the FIREFLy (orange histograms) and Juniper (green histograms) data reductions. The median (solid lines) and 1σ credible intervals for each parameter (dashed lines) are overlaid. c, Residuals between the data and best-fitting model.

Besides hydrocarbons, we detect thermal emission from the observer-facing nightside hemisphere ($18.5\sigma /\mathrm{ln}{\mathcal{B}}=169$ for Juniper and $17.3\sigma /\mathrm{ln}{\mathcal{B}}=146$ for FIREFLy) and aerosols ($5.7\sigma /\mathrm{ln}{\mathcal{B}}=14.5$ for Juniper and $5.3\sigma /\mathrm{ln}{\mathcal{B}}=12.3$ for FIREFLy). Emission from the nightside hemisphere (about 400 K) causes the substantial downwards slope longwards of 3.5 μm. Because thermal emission in opacity windows arises from deeper, hotter layers, an optically thick cloud deck must exist for pressures deeper than roughly 10 mbar (Fig. 2c) to explain the lack of ‘negative spikes’ between the CH₄ absorption bands (Fig. 3a,c). Because most thermal emission arises from the cloud-top temperature, the recent measurement of an approximately 190 K thermal excess in MIRI photometry of the WD 1856 system¹² (at a different epoch near quarter phase) potentially indicates variable cloud-top pressures for different hemispheres or in time. Alternatively, an unidentified systematic or calibration issue affecting either instrument may reconcile the difference, but we did not identify any such effect in either dataset. Further evidence of aerosols arises from the scattering slope shortwards of 1 μm, which is more prominent than H₂ Rayleigh scattering from a clear atmosphere. We investigated several specific aerosols using a Mie scattering model, but the present data are insufficient to identify the chemical composition of the aerosol (Methods). Finally, we report the first constrained measurement of the mass of WD 1856 b: 4.3–10.9 M_J (1σ spread across both data reductions or $6.{7}_{-2.4}^{+2.8}\,{M}_{{\rm{J}}}$ for FIREFLy and $7.{8}_{-2.7}^{+3.1}\,{M}_{{\rm{J}}}$ for Juniper). Our mass constraint arises from a scale height trade between the CH₄ absorption features and the integrated optical depth to the emitting pressure of the nightside atmosphere.

The atmospheric temperature at which WD 1856 b radiates to space (T_eff = 390–412 K, or $40{0}_{-10}^{+6}\,{\rm{K}}$ for FIREFLy and $40{5}_{-11}^{+7}\,{\rm{K}}$ for Juniper; Methods) is substantially increased over both the equilibrium temperature (160 K) and that expected for an evolved giant planet at the roughly 10 Gyr system age (≲100 K)¹. Given the retrieved mass and the present-day circular orbit, internal power sources such as deuterium fusion or tidal heating cannot contribute to the observed effective temperature (Methods). However, a thermal reset of WD 1856 b can have occurred by either tidal heating during high-eccentricity migration¹⁰ or immersion in the stellar envelope during a common-envelope phase^8,9. These two scenarios predict different migration times relative to the death of the host star: common-envelope evolution coincides with the end of the AGB phase of the progenitor (5.4 ± 0.7 Gyr ago; Methods) and lasts about 1 Myr (ref. ²⁰), whereas high-eccentricity migration can occur any time during the white dwarf phase. Because substellar objects cool down at a predictable rate, our planetary mass and temperature constraints provide critical information to infer the reheating epoch.

We reconstructed the thermal evolution of WD 1856 b using theoretical cooling models for substellar objects (Methods) to extrapolate its effective temperature backwards from the present. Figure 4 shows a random set of reconstructed thermal histories for WD 1856 b, alongside the equilibrium temperature evolution (tracing the cooling of the white dwarf). The range of possible histories is governed by uncertainties on the mass of WD 1856 b and the cooling age of the host. For each thermal history, we calculated the cooling age of the white dwarf at the time of reheating (denoted t₀; Methods). In Fig. 4, t₀ corresponds roughly to the time when each cooling model intercepts the top of the diagram. We find that the reheating of WD 1856 b occurred 3.0–5.5 Gyr after the end of the AGB phase ($4.{2}_{-1.2}^{+1.0}\,{\rm{Gyr}}$ for FIREFLy and $4.{6}_{-1.0}^{+0.9}\,{\rm{Gyr}}$ for Juniper; Fig. 4b). More conservative 2σ lower bounds on t₀ imply that reheating occurred at least 1.4 Gyr (FIREFLy) or 2.1 Gyr (Juniper) after the AGB phase. The AGB and post-AGB/pre-white-dwarf phases are extremely brief by comparison, lasting less than 2 Myr and 0.1 Myr, respectively (Methods). The timing of the inferred reheating event is, therefore, inconsistent with common-envelope evolution during either phase. Therefore, WD 1856 b most likely underwent high-eccentricity migration to its current orbit, with the inferred reheating event corresponding to tidal circularization.

Fig. 4: Evolutionary history of WD 1856 b.

a, Reconstructed thermal evolution of WD 1856 b. An ensemble of randomly drawn thermal histories compatible with the JWST-derived constraints on the mass and effective temperature of WD 1856 b and the cooling age of the host white dwarf are shown (orange curves for the FIREFLy data, green curves for Juniper). The thermal history for an object at the 0.02-au orbit of WD 1856 b maintaining the zero-albedo equilibrium temperature is overlaid for comparison (dashed black curve). b, Inferred host cooling age at the time of reheating of WD 1856 b during migration, t₀, based on the ensemble of reconstructed thermal histories (orange and green histograms for FIREFLy and Juniper, respectively). The median (solid line) and 1σ credible interval (dashed lines) are overlaid and also plotted above panel a. Because common-envelope evolution predicts t₀ consistent with zero, the present-day mass and temperature favour migration billions of years after the conclusion of the AGB phase.

WD 1856 b represents the first well-characterized transiting planet orbiting a white dwarf. The inferred thermal evolution of WD 1856 b demonstrates that high-eccentricity migration is a plausible fate for giant planets after the stellar main sequence. The retrieved CH₄ abundance is similar to the deep atmosphere of Neptune (4% (ref. ²¹)), which requires notable carbon enrichment of the planet’s H₂ envelope from volatile-rich material, whether accreted before its migration^22,23,24 or afterwards²⁵. This high atmospheric metallicity (≳100 × solar) enhances aerosol production^26,27,28, consistent with the detection of a short-wavelength scattering slope in our transmission spectrum. As WD 1856 b demonstrates, spectroscopy of planets orbiting white dwarfs offers a new opportunity to determine the fate of planetary systems after the death of their star.

Methods

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:

$${A}_{{\rm{p}}}(d)={R}_{{\rm{p}}}^{2}\theta +{R}_{\ast }^{2}\phi -\frac{1}{2}{R}_{{\rm{p}}}^{2}\sin (2\theta )-\frac{1}{2}{R}_{\ast }^{2}\sin (2\phi )$$

(1)

in which

$$\theta ={\cos }^{-1}\,\left(\frac{{d}^{2}+{R}_{{\rm{p}}}^{2}-{R}_{\ast }^{2}}{2d{R}_{{\rm{p}}}}\right)$$

(2)

$$\phi ={\cos }^{-1}\,\left(\frac{{d}^{2}+{R}_{\ast }^{2}-{R}_{{\rm{p}}}^{2}}{2d{R}_{\ast }}\right)$$

(3)

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:

$$\frac{{A}_{{\rm{p}}}}{{A}_{\ast }}=\frac{1}{{\rm{\pi }}}\left[{\left(\frac{{R}_{{\rm{p}}}}{{R}_{\ast }}\right)}^{2}\left(\theta -\frac{1}{2}\sin 2\theta \right)+\left(\phi -\frac{1}{2}\sin 2\phi \right)\right]$$

(4)

in which

$$\theta ={\cos }^{-1}\,\left[\frac{{b}^{2}+{\left(\frac{{R}_{{\rm{p}}}}{{R}_{\ast }}\right)}^{2}-1}{2b\left(\frac{{R}_{{\rm{p}}}}{{R}_{\ast }}\right)}\right]$$

(5)

$$\phi ={\cos }^{-1}\,\left[\frac{{b}^{2}-{\left(\frac{{R}_{{\rm{p}}}}{{R}_{\ast }}\right)}^{2}+1}{2b}\right]$$

(6)

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³⁹:

$${\varDelta }_{\lambda }=\frac{{A}_{{\rm{p}}({\rm{t}}{\rm{o}}{\rm{p}})}-{\int }_{{A}_{{\rm{p}}}}{{\mathcal{T}}}_{\lambda }{\rm{d}}A}{{\rm{\pi }}{R}_{\ast }^{2}}\left(\frac{1}{1+\frac{{F}_{{\rm{p}}({\rm{n}}{\rm{i}}{\rm{g}}{\rm{h}}{\rm{t}}),\lambda }}{{F}_{\ast ,\lambda }}}\right)$$

(7)

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:

$${\varDelta }_{\lambda }=\frac{{A}_{{\rm{p}}({\rm{t}}{\rm{o}}{\rm{p}})}-{\int }_{{A}_{{\rm{p}}}}{{\mathcal{T}}}_{\lambda }{\rm{d}}A}{{\rm{\pi }}{R}_{\ast }^{2}}\times \left(\frac{1}{1+\frac{{R}_{{\rm{p}},({\rm{n}}{\rm{i}}{\rm{g}}{\rm{h}}{\rm{t}}),\lambda }^{2}}{{R}_{\ast }^{2}}\frac{{F}_{{\rm{p}}({\rm{n}}{\rm{i}}{\rm{g}}{\rm{h}}{\rm{t}}),{\rm{s}}{\rm{u}}{\rm{r}}{\rm{f}},\lambda }}{{F}_{\ast ,{\rm{s}}{\rm{u}}{\rm{r}}{\rm{f}},\lambda }}}\right)$$

(8)

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:

$${\mathcal{R}}=4{\left(\frac{{T}_{{\rm{eff}},{\rm{p}}}}{{T}_{{\rm{eff}},\ast }}\right)}^{4}{\left(\frac{{a}_{{\rm{orb}}}}{{R}_{\ast }}\right)}^{2},$$

(9)

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:

$$\begin{array}{c}{T}_{{\rm{e}}{\rm{f}}{\rm{f}},{\rm{p}}}={\left(\frac{21}{8{\rm{\pi }}}\frac{{G}^{2}{M}_{\ast }^{3}{R}_{{\rm{p}}}^{3}}{{\sigma }_{{\rm{S}}{\rm{B}}}{a}_{{\rm{o}}{\rm{r}}{\rm{b}}}^{9}}{k}_{2{\rm{p}}}{\tau }_{{\rm{p}}}{e}_{{\rm{o}}{\rm{r}}{\rm{b}}}^{2}\right)}^{1/4}\\ \,\approx \,400\,{\rm{K}}{\left(\frac{{M}_{\ast }}{0.60{M}_{\odot }}\frac{{R}_{{\rm{p}}}}{0.91{R}_{{\rm{J}}}}\right)}^{3/4}{\left(\frac{{a}_{{\rm{o}}{\rm{r}}{\rm{b}}}}{0.02{\rm{A}}{\rm{U}}}\right)}^{-9/4}\\ \,\times \,{\left(\frac{{k}_{2{\rm{p}}}}{0.4}\frac{{\tau }_{{\rm{p}}}}{0.1{\rm{s}}}\right)}^{1/4}{\left(\frac{{e}_{{\rm{o}}{\rm{r}}{\rm{b}}}}{0.02}\right)}^{1/2}.\end{array}$$

(10)

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.