SiO and a super-stellar C/O ratio in the atmosphere of the giant exoplanet WASP-121 b

Abstract

Refractory elements such as iron, magnesium and silicon can be detected in the atmospheres of ultrahot giant planets. This provides an opportunity to quantify the amount of refractory material accreted during formation, along with volatile gases and ices. However, simultaneous detections of refractories and volatiles have proved challenging, as the most prominent spectral features of associated atoms and molecules span a broad wavelength range. Here, using a single JWST observation of the ultrahot giant planet WASP-121 b, we report detections of H₂O (5.5–13.5σ), CO (10.8–12.8σ) and SiO (5.7–6.2σ) in the planet’s dayside atmosphere and CH₄ (3.1–5.1σ) in the nightside atmosphere. We measure super-stellar values for the atmospheric C/H, O/H, Si/H and C/O ratios, which point to the joint importance of pebbles and planetesimals in giant planet formation. The CH₄-rich nightside composition is also indicative of dynamical processes, such as strong vertical mixing, having a profound influence on the chemistry of ultrahot giant planets.

Main

WASP-121 b is an ‘ultrahot’ giant planet with dayside temperatures exceeding ~2,000 K, high enough to prevent refractory elements from condensing. Provided the refractories are not cold-trapped on the cooler nightside, this makes it possible to constrain the refractory and volatile contents of the atmospheric envelope by measuring the dayside gas composition¹. When combined with key elemental ratios such as the C/O ratio, the derived refractory and volatile abundances can provide insights into the planet’s formation history². A practical challenge is that the wavelength coverages of available spectrographs have typically allowed detections of either refractory or volatile species, but not both, with the same observation. This is because the strongest spectral features of major refractory species—including Fe, Mg, Ca, SiO, TiO and VO—are at near-ultraviolet (near-UV) and optical wavelengths, whereas those of major volatile species such as H₂O and CO occur at infrared wavelengths^3,4.

We have used JWST to observe a ‘phase curve’ for WASP-121 b, in which thermal emission from the system was continuously recorded as the planet completed an orbit around the host star⁵. As expected, the infrared wavelength range spanned by the measurement provides sensitivity to important volatiles such as H₂O and CO in the planetary atmosphere. The precision of JWST also allows us to detect a relatively weak SiO band, delivering abundance constraints for a key refractory species. In addition, the complete orbital phase coverage enables separate characterization of the planet’s dayside and nightside temperature profiles and chemical compositions. This, in turn, allows us to explore how dynamical processes shape the global atmospheric chemistry, which can be driven out of equilibrium with local thermobaric conditions if gas transport occurs on sufficiently short timescales^6,7,8,9.

Results

Observations

We observed WASP-121 b using the Near Infrared Spectrograph (NIRSpec) instrument on the James Webb Space Telescope (JWST) with the G395H grism, covering the wavelength (λ) range 2.73–5.17 μm. The observation lasted 37.8 h and encompassed two consecutive secondary eclipses and a primary transit. We generated 349 wavelength-dependent ‘spectroscopic’ phase curves by binning each measured spectrum into 10-pixel-wide bins along the dispersion axis, translating to spectroscopic bin widths of Δλ = 6.7 nm and a mean spectral resolving power of λ/Δλ ≈ 600 across the G395H passband. Further details of the observation and data reduction are provided in ‘Observations and data reduction’ in Methods.

Phase-curve fitting

To fit the spectroscopic phase curves, we adopted a spherical harmonic dipole to describe the planetary brightness map and a quadratic limb-darkening law for the host star. A linear polynomial in time and the pixel coordinates of the target spectrum on the detector was simultaneously fitted to account for drift in the instrument baseline. We extracted the following three observables to characterize the planetary atmosphere from each spectroscopic phase-curve fit: the ratio of the disk-integrated planetary brightness to the disk-integrated stellar brightness (F_p/F_⋆) as a function of orbital phase (ɸ); the phase offset of the planetary brightness map (Δɸ); and the effective planetary radius (R_p). Further details of the phase-curve fitting are given in ‘White phase-curve fitting’ and ‘Spectroscopic phase-curve fitting’ in Methods.

We binned the measured F_p/F_⋆ values into 36 phase bins that avoided the two eclipses and transit, with each phase bin covering 46 min of the planetary orbit (Extended Data Fig. 1). By doing this for each spectroscopic phase curve, we were able to obtain the planetary emission spectrum as a function of orbital phase (see ‘Generating the phase-resolved emission spectra’ in Methods). Figure 1 shows the resulting emission spectra for the phase bins immediately preceding the secondary eclipse and primary transit, when the dayside and nightside hemispheres of WASP-121 b were maximally visible, respectively. To verify the robustness of these spectra, we confirmed that they were in good agreement with those produced by a second independent analysis (see ‘Independent data reduction’ in Methods). Emission spectra for the full set of orbital phase bins are shown in Extended Data Fig. 2.

Fig. 1: Dayside and nightside emission spectra for WASP-121 b.

a, Dayside planet-to-star emission based on n = 5,000 posterior samples from each of the 349 spectroscopic phase-curve fits. The grey circles show median values and the grey error bars indicate 1σ uncertainties. The black circles show the median F_p/F_⋆ values in 50 equally spaced bins, with black error bars showing ${\widetilde{\sigma }}_{i}/\sqrt{{N}_{i}}$, where ${\widetilde{\sigma }}_{i}$ is the median 1σ uncertainty and N_i is the number of data points in the ith bin. The orange line shows the maximum a posteriori model that assumes thermochemical equilibrium, with detected molecules labelled below. Predicted dayside spectra from the GCM simulations published in ref. ¹¹ for atmospheric metallicities of 1 × solar (light blue line) and 5 × solar (dark blue line) are also shown. b, The same as a for the nightside spectrum. The H₂O and CO absorption features in the GCM spectra are labelled but not seen in the data. The dashed red line shows the spectrum predicted by the maximum a posteriori model (orange line) when CO opacity is switched off. c, Wavelength-dependent brightness temperatures (T_b) derived from the binned dayside emission data and model shown in a. d, Wavelength-dependent T_b derived from the binned nightside emission data and models shown in b.

Atmospheric retrieval analyses

To infer the atmospheric properties of the dayside and nightside hemispheres, we performed retrieval analyses for the two spectra shown in Fig. 1. We adopted widely used parametric forms for the dayside and nightside pressure–temperature (PT) profiles and assumed thermochemical equilibrium. We allowed the elemental abundances of carbon, oxygen and silicon to vary individually as free parameters, following tests that indicated molecules comprising these atoms accounted for the spectral features in the data. We set the relative abundances of all other heavy elements equal to the solar values and allowed them to vary collectively in the retrievals via an additional metallicity parameter. Further details are given in ‘Atmospheric retrieval analyses’ in Methods.

The outputs of these retrieval analyses are shown in Fig. 2. The PT profile for the dayside hemisphere exhibits a thermal inversion, in line with previous measurements for WASP-121 b¹⁰. We find that the temperature increases from approximately 2,600 to 3,200 K over the 10⁻¹ to 10⁻³ bar pressure range. Across the same pressure range, the inferred nightside PT profile cools from approximately 1,400 to 1,100 K. To interpret the retrieved abundances for carbon, oxygen and silicon, we compared them to updated elemental abundances of the host star obtained from high-resolution spectroscopy (see ‘Stellar abundances and relative planetary abundances’ in Methods). For the dayside spectrum, we measure planetary ratios relative to the equivalent host-star ratios of ${({\rm{C}}/{\rm{H}})/({\rm{C}}/{\rm{H}})}_{\star }={23.96}_{-4.75}^{+6.13}$, ${({\rm{O}}/{\rm{H}})/({\rm{O}}/{\rm{H}})}_{\star }={12.19}_{-2.24}^{+2.78}$ and ${({\rm{Si}}/{\rm{H}})/({\rm{Si}}/{\rm{H}})}_{\star }={9.89}_{-2.89}^{+5.99}$, where quoted values are the medians with 1σ uncertainties defined by the 16th and 84th percentiles (68% credible interval). The abundances for all three elements are super-stellar, with (C/H)/(C/H)_⋆ > 11.92, (O/H)/(O/H)_⋆ > 7.27 and (Si/H)/(Si/H)_⋆ > 2.49 at 99.99% probability. We derive a C/O ratio of ${0.92}_{-0.03}^{+0.02}$, which is >1.63 times the stellar ratio of (C/O)_⋆ = 0.47 ± 0.06 at 99.99% probability. We also determine upper limits for the volatile/refractory ratios of (C/Si)/(C/Si)_⋆ < 8.07, (O/Si)/(O/Si)_⋆ < 4.02 and ((C + O)/Si)/((C + O)/Si)_⋆ < 5.24 at 99.99% probability. Posterior distributions for these elemental abundances and the derived ratios are shown in Fig. 3.

Fig. 2: Properties of the atmosphere of WASP-121 b.

a, PT profiles inferred from the retrieval analysis for the dayside (red) and nightside (blue) atmospheres assuming thermochemical equilibrium. Condensation curves for refractory species are included as labelled dash-dotted lines. b, Normalized contribution functions for the dayside (red) and nightside (blue) atmospheres integrated across the G395H passband, illustrating the relative amount of radiation originating from each pressure level. c, Volume-mixing ratios (VMRs) inferred for key molecules in the nightside atmosphere. d, The same as c for the dayside atmosphere. Solid lines in a, c and d are the median values obtained at each pressure level and the shaded regions indicate the associated 1σ uncertainties defined by the 16th and 84th percentiles (68% credible intervals) based on n = 1,000 posterior samples.

Fig. 3: Elemental abundances and derived ratios.

Top: posterior distributions obtained for C/H (left), O/H (middle) and Si/H (right) ratios from the retrieval analyses for the dayside atmosphere assuming thermochemical equilibrium. Bottom: posterior histograms derived for the C/O ratio (left), volatile/refractory ratios (middle) and metallicity (right). The metallicity regulated the abundances of heavy elements besides C, O and Si, as described in the text. Vertical dashed lines indicate the corresponding values measured for the host star using high-resolution spectroscopy.

The abundances for different chemical species implied by the retrievals are shown in Fig. 2 and the corresponding contributions to the dayside and nightside emission spectra are illustrated in Extended Data Fig. 3. On the dayside, the data trace out spectral features of H₂O, SiO and CO in emission, with the H₂O and SiO features muted considerably owing to thermal dissociation (Fig. 2). The measured amplitudes of the H₂O and CO dayside emission features are comparable to those predicted by three-dimensional (3D) general circulation model (GCM) simulations for WASP-121 b published in ref. ¹¹ (Fig. 1), although the GCM spectra do not show an SiO feature. On the nightside, the data exhibit spectral features of CH₄ in absorption, rather than H₂O in absorption as had previously been reported for WASP-121 b on the basis of HST observations at shorter wavelengths¹². Using ‘free chemistry’ retrieval analyses that did not enforce the requirement of thermochemical equilibrium (see ‘Atmospheric retrieval analyses’ in Methods), we obtain statistically significant detections of dayside H₂O (5.5–13.5σ), dayside CO (10.8–12.8σ), dayside SiO (5.7–6.2σ) and nightside CH₄ (3.1–5.1σ).

Discussion

We have obtained a conclusive detection of SiO in the atmosphere of WASP-121 b. Evidence of SiO had previously been reported in the atmosphere of WASP-121 b and another ultrahot gas giant, WASP-178 b, on the basis of absorption signals measured at near-UV wavelengths¹³. However, these near-UV data could also be explained by Mg i and Fe ii absorption, preventing an unambiguous identification of SiO. With NIRSpec, we are now able to verify the presence of SiO in the atmosphere of WASP-121 b and, due to its properties as a strong absorber of shortwave radiation¹³, the important role that it is likely playing in driving the dayside thermal inversion. Furthermore, silicon is expected to be one of the primary cloud constituents in hot (~1,000–2,000 K) atmospheres¹⁴ and can provide valuable information about the overall formation history of a planet, as it allows the volatile/refractory ratio to be constrained^1,2. By delivering a robust measurement of a planetary atmosphere’s silicon abundance, our SiO detection stands out from previous detections of atomic silicon^15,16 and silicate clouds^17,18 that lacked abundance constraints.

A GCM simulation for WASP-121 b that assumed solar metallicity and included SiO opacity was presented in ref. ¹⁹. The predicted dayside spectrum did not show the 4–4.3 μm SiO band, as it was obscured by overlapping H₂O and CO bands. Our NIRSpec observation demonstrates that in reality this will not always be the case, particularly if an atmosphere is enriched in silicon and oxygen, as we find for WASP-121 b. It raises the possibility of further SiO detections with JWST in the future, which would complement efforts at UV wavelengths while being far less hampered by interstellar dust extinction and low UV fluxes of host stars.

Our finding that the C/H, O/H and C/O ratios of WASP-121 b are all super-stellar provides evidence that the planet accreted most of its volatiles as gas that was enriched by inward-drifting pebbles^20,21,22. The super-stellar C/O and C/H ratios in particular may indicate that WASP-121 b accreted most of its gaseous envelope beyond the H₂O ice line but within the CH₄ ice line, where the gas-phase carbon abundance could have been locally elevated by the evaporation of inward-drifting CH₄-rich pebbles²². If WASP-121 b instead accreted most of its volatiles as icy planetesimals between the H₂O and CO ice lines, the oxygen-rich composition of these planetesimals should have produced a sub-stellar C/O ratio²³. Alternatively, a sub-stellar volatile inventory would be expected if the volatiles were accreted from gas that had not been enriched by evaporating pebbles²¹.

While gas accretion seems to have delivered the bulk of volatiles, our finding of a super-stellar Si/H ratio (Fig. 3) suggests that there was also substantial enrichment by refractory solids following accretion of the gaseous envelope^24,25. We calculate that our measured silicon abundance is consistent with the addition of approximately ${21.2}_{-6.6}^{+11.8}\,M_{\oplus}$ (where $M_{\oplus}$ is the mass of the Earth) of rocky material to the atmosphere (see ‘Accretion of rocky material’ in Methods), which is comparable to predictions of planetesimal accretion models^26,27. By contrast, it may not be possible for pebbles to deliver this quantity of refractory material, as the accretion of solid pebbles should largely cease once the pebble isolation mass is reached²⁸. After this point, the planet would continue to accrete volatile-rich gas and obtain an atmosphere that has a sub-stellar silicon abundance or a higher volatile/refractory ratio than we have observed (Fig. 3). For example, in the pebble accretion simulations of ref. ²⁴ that produced super-stellar C/H, O/H and C/O values, the (C + O)/Si ratio rises to at least ten times the stellar ratio.

Super-stellar C/O constraints have recently been obtained for WASP-121 b using ground-based spectroscopy^3,4, underscoring the robustness of this particular result. Combined with their own measurement of a super-stellar volatile/refractory ratio, the authors of ref. ⁴ proposed a broadly similar formation scenario to the one that we have outlined, with WASP-121 b accreting most of its atmospheric envelope in a volatile-rich environment beyond the H₂O ice line. In contrast, the authors of ref. ³ recovered a sub-stellar volatile/refractory ratio, which they interpreted as evidence for formation within the H₂O ice line and accretion dominated by refractory-rich planetesimals rather than volatile ices. Ultimately, these tensions between available datasets will need to be resolved to uncover a clearer picture of how WASP-121 b formed.

We also stress the broader difficulties in linking a planet’s observed composition to its formation history. Various complications are reviewed in refs. ^1,2,22,29 and include: accretion from compositionally distinct regions within the protoplanetary disk due to orbital migration; unknown dust opacity and self-shadowing affecting the thermal structure and associated chemistry of the disk; time-dependent chemical and physical evolution of the disk; and erosion of the planetary core raising the heavy element content of the atmospheric envelope. Further investigation of such effects will be required to rigorously assess the basic qualitative picture we have presented here for the formation of WASP-121 b.

Regarding the chemistry of the atmosphere, the evidence we have uncovered for CH₄ on the nightside of WASP-121 b was not anticipated. One reason for this is that models of hot giant planet atmospheres have tended to assume atmospheric C/O ratios close to the solar value of approximately 0.55 (ref. ³⁰). This leads to dayside and nightside spectra sculpted by H₂O and CO, as seen for the GCM of ref. ¹¹ (Fig. 1). When the C/O ratio is increased from 0.55 to our measured value of 0.92, the nightside CO abundance changes little, but the CH₄ abundance increases by 2–4 orders of magnitude and the H₂O abundance decreases by a similarly large amount (Extended Data Fig. 4). Sequestration of oxygen into the various oxygen-bearing condensates shown in Fig. 2 can further reduce the H₂O abundance³¹. As a result, the atmospheric opacity becomes dominated by CH₄, explaining why H₂O, CO and SiO are not detected on the nightside (Extended Data Fig. 3). Although the CO opacity is stronger than the CH₄ opacity within narrow line cores across the 4–4.5 μm wavelength range, the effect this has on the emission spectrum is too subtle to be detected at the resolution of the data (Fig. 1).

When the dynamics of the atmosphere are taken into account, a high C/O ratio by itself cannot easily explain our measured abundance of CH₄ on the nightside of WASP-121 b. Simulations of hot giant planets find that horizontal winds with speeds of ~1–10 km s⁻¹ advect gas between the dayside and nightside hemispheres at pressures below ~1 bar (refs. ^8,19,32,33). This would mean that gas traverses the nightside hemisphere of WASP-121 b on typical timescales of ~10⁴–10⁶ s, which is substantially shorter than the timescales required to reach chemical equilibrium. Specifically, we find that the nightside infrared photosphere coincides with pressures of ~10²–10⁻⁴ bar and temperatures of ~1,200 K (Fig. 2), where the timescales for CH₄ and CO to reach their local equilibrium abundances would be ~10¹⁰ s and ~10¹⁵ s, respectively^9,34. Conversely, on the dayside, high temperatures translate to chemical reaction timescales that are shorter than dynamical timescales throughout the observable atmosphere, allowing chemical equilibrium to be established. If horizontal winds dominate the transport of gas throughout the observable atmosphere⁸, we would expect gas to reach equilibrium with the dayside conditions and then be quenched as it travels across the nightside, preserving its CH₄-poor composition.

We propose instead that vertical mixing on the nightside hemisphere is efficiently transporting CH₄-rich gas up to the infrared photosphere from deeper layers of the atmosphere that are in chemical equilibrium. To test this hypothesis, we performed 1D model calculations for the nightside chemistry with vertical mixing parameterized by an eddy diffusion coefficient (K_zz). We adopted the retrieved nightside PT profile shown in Fig. 2 but started with the dayside chemical composition. The latter can be thought of as a conservative assumption, equivalent to instantaneous advection of gas from the dayside to the nightside. We then evolved the composition of the gas over timescales ranging from 10⁴–10⁶ s for a range of K_zz values between 10⁹ cm² s⁻¹ and 10¹³ cm² s⁻¹ using a C–H–N–O–S non-equilibrium chemical network³⁵. Figure 4 shows the resulting CH₄ abundances predicted at the infrared photosphere. We find that the measured CH₄ abundance can be reproduced with K_zz values of ~10⁹–10¹² cm² s⁻¹, which are consistent with available constraints^36,37. The predicted CH₄ abundance is also higher than the H₂O abundance, in line with our detection of CH₄ and non-detection of H₂O on the nightside. Similar results were obtained when we repeated the calculations using the PT profiles from the GCM simulation of ref. ¹¹ (Extended Data Fig. 5). From these analyses, we deduce that gas mixed vertically from deeper layers of the atmosphere is setting the composition at the infrared photosphere, rather than horizontally quenched CH₄-poor gas transported from the dayside, as is often seen in dynamical simulations^7,8.

Fig. 4: Effect of vertical mixing on nightside abundances.

a, Pressure-dependent VMR profiles for CH₄, H₂O and CO assuming equilibrium chemistry (dashed lines) and non-equilibrium chemistry (solid lines). The non-equilibrium chemistry calculations adopted K_zz = 10¹¹ cm² s⁻¹ and were evolved over a timescale of 10⁵ s. The light grey shading indicates the pressure levels probed by the data. Dark orange shading shows the range of median CH₄ VMRs obtained from the retrieval analyses and light orange shading shows the combined extent of the associated 1σ uncertainties defined by the 16th and 84th percentiles. b, Mean CH₄ VMR across the 10⁻¹–10⁻³ bar pressure range as a function of K_zz over the timescales indicated. The orange shading is the same as in a. The dashed grey line shows the retrieved CH₄ VMR for the dayside, which was adopted as the starting abundance.

Vertical mixing also offers a natural explanation for why silicon is not cold-trapped on the nightside hemisphere of WASP-121 b, as demonstrated by our detection of SiO on the dayside. Recent microphysical modelling has found that cloud formation should occur rapidly in hot giant planet atmospheres on timescales of ~1 s (ref. ³⁸). We thus expect clouds of silicates to form efficiently on the nightside of WASP-121 b at pressures of ~10⁻¹–1 bar, where the measured temperature profile crosses the condensation curves of MgSiO₃, Mg₂SiO₄ and SiO₂ (Fig. 2). As mentioned above, iron and magnesium have also been detected in the atmosphere of WASP-121 b, along with numerous other refractory metals such as vanadium, chromium and nickel^39,40,41. The same mixing that we hypothesize is transporting CH₄ upwards to the observable atmosphere is likely preventing these refractory species from condensing on the nightside and becoming permanently cold-trapped in the deep atmosphere. A notable exception is titanium, for which non-detections provide strong evidence of being efficiently cold-trapped^40,41. However, despite conditions being favourable for the formation of various refractory clouds on the nightside, we find that models including clouds do not appreciably improve the fit to the available data (see ‘Discussion of the retrieval analyses’ in Methods). At the same time, we are unable to rule out clouds at pressures above ~10⁻¹ bar, as such clouds would not necessarily produce an observable signature (Extended Data Fig. 6).

GCM simulations will be valuable for further investigating the potential role played by nightside clouds. For example, existing cloud-free GCMs predict nightside temperatures that are warmer than we measure, particularly at pressures >10⁻² bar (Extended Data Fig. 5), resulting in strong emission peaks between absorption bands that are not observed in the data (Fig. 1). Nightside clouds may help to reduce this discrepancy by blocking emission from the deepest, warmest layers of the atmosphere^42,43,44,45. Other effects, such as magnetic influences through Lorentz forces, could also play an important role in lowering nightside temperatures by reducing the efficiency of day–night heat recirculation^46,47,48.

Methods

Observations and data reduction

The WASP-121 system was continuously monitored for 37.8 h on 2022 October 14–15 using the Bright Object Time Series mode of the NIRSpec instrument on JWST. The observation began 145 min before secondary eclipse ingress and finished 105 min after egress of the following secondary eclipse. As the planet is obscured during secondary eclipses, this allowed us to calibrate the brightness level of the host star at both the start and end of the observation. Light received from the WASP-121 b system was dispersed across the two 2,048 × 2,048 pixel HgCdTe detectors at the NIRSpec focal plane by the G395H grism. The two detectors, denoted NRS1 and NRS2, provide continuous wavelength coverage between 2.70–3.72 μm and 3.82–5.15 μm, respectively. For each integration, 2,048 × 32 pixel subarrays were read from both detectors using the NRSRAPID readout pattern with 42 groups per integration, translating to individual integration times of 38.8 s. A total of 3,504 integrations was acquired over the full observation. There were five high-gain antenna movements during the observation, which caused small, but noticeable, deviations in the measured fluxes for five of the integrations. We retained these impacted integrations in our analysis, as we expect them to have a negligible effect on the overall results.

To analyse the data, we started by using the FIREFLy code⁴⁹ to extract the source spectrum from each integration. Further details of our FIREFLy reduction can be found in ref. ⁵. From the resulting time series of spectra, separate ‘white’ phase curves were generated for the NRS1 and NRS2 detectors by summing the spectra from each detector across the full wavelength range. We also binned the spectra into narrower 10-pixel-wide channels along the dispersion axis, producing 349 spectroscopic phase curves.

White phase-curve fitting

Before fitting the spectroscopic phase curves, we refined our estimates for various system parameters by performing a joint fit to the higher-signal-to-noise white phase curves. For each white phase curve, our model Π took the form Π = α × β, where α is the astrophysical signal and β is the instrument baseline.

We used the starry package⁵⁰ to model the time-dependent astrophysical signal α. starry accounts for emission from both the star and the planet, as well as primary transits and secondary eclipses. We adopted a degree ${\mathcal{l}}=1$ spherical harmonic expansion for the planetary brightness map and a quadratic limb-darkening profile for the star.

For the instrument baseline β, we investigated three forms. First, a linear trend with respect to time t:

$${\beta}_{1}\left(t\right)={c}_{0}+{c}_{1}t$$

Second, a linear trend with respect to t and the xy pointing jitter coordinates described in ref. ⁵:

$${\beta }_{2}\left(t,x,y\right)={c}_{0}+{c}_{1}t+{c}_{x}x+{c}_{y}y.$$

Third, a linear trend with respect to t, plus an exponential function:

$${\beta }_{3}\left(t\right)={c}_{0}+{c}_{1}t+{r}_{0}{\mathrm{e}}^{-t/{r}_{1}}.$$

We note that β₃ allows for an initial settling of the telescope following a new pointing. Baselines of this form have mostly been used for analyses of JWST Mid-Infrared Instrument (MIRI) light curves, which tend to exhibit prominent ramp-like systematics lasting ~30–60 min at the beginning of an exposure⁵¹.

In the fitting, the following parameters were shared across both the NRS1 and NRS2 white phase curves: stellar mass (M_⋆); planet mass (M_p); orbital inclination (i); logarithm of the stellar radius (log₁₀R_⋆); and deviation in the orbital period (ΔP_orb) from the value of P_orb = 1.27492504 ± 1.5 × 10⁻⁷ days reported in ref. ⁵². We fitted the following parameters separately for each detector: the logarithm of the planetary radius (log₁₀R_p); the deviation in the transit mid-time (ΔT_mid) from the value of T_mid = 2,459,867.64297 ± 0.00027 BJD_TDB predicted by the ephemerides of ref. ⁵²; the relative planet-to-star luminosity (A); the coefficient of the Y_1,0 spherical harmonic (y_1,0); the phase offset of the planetary brightness map (Δɸ); the coefficients of the stellar limb-darkening profile (u₁ and u₂); the instrument baseline parameters (c₀, c₁, c_x, c_y, r₀ and r₁); and a high-frequency systematics noise term (σ_syst), explained below. We also tested more elaborate planetary brightness maps comprising combinations of the Y_1,0, Y_1,1 and Y_2,0 spherical harmonics, with the associated coefficients treated as free parameters in the model fitting. However, we found that a planetary brightness map described by Y_1,0 alone was favoured by the Bayesian information criterion (BIC).

For each white phase-curve fit, a multivariate normal distribution was adopted for the likelihood function and marginalized using affine-invariant Markov chain Monte Carlo, as implemented by the emcee package⁵³. We used 150 walkers and took 5,000 steps to obtain samples from the posterior distribution, following an initial burn-in phase of 500 steps. Note that when evaluating the likelihood function, the error bar on the ith data point was taken to be the quadrature sum ${\sigma }_{i}=\sqrt{{\sigma }_{\text{Poiss},i}^{2}+{\sigma }_{\text{syst}}^{2}}$, where σ_Poiss,i was the estimated Poisson noise for the ith data point returned by the instrument pipeline and σ_syst was a free parameter allowing for additional high-frequency systematics.

This overall methodology for fitting the white phase curves is similar to that used in ref. ⁵, but with a few variations. First, we adopted different priors for a subset of parameters. In particular, normal priors were adopted for M_⋆ and M_p based on a radial velocity analysis of the G395H data⁵⁴. Normal priors obtained using the exoTiC-LD package⁵⁵ were also adopted for the wavelength-dependent limb-darkening coefficients (u₁ and u₂). Details of all adopted priors are provided in Supplementary Fig. 1. Second, we allowed ΔT_mid to vary separately for the NRS1 and NRS2 detectors, whereas ref. ⁵ adopted a shared transit mid-time for both detectors. Third, although the beginning of the exposure was not visibly affected by strong systematics, we discarded the first 45 min of data before fitting. This precautionary step was taken to provide ample time for the telescope to settle into its new pointing, while still retaining 100 min of pre-eclipse baseline.

Of the three instrumental baselines tested (β₁, β₂ and β₃), we found that β₂ gave the lowest BIC. The results of the corresponding fit are summarized in Supplementary Table 1 and Supplementary Fig. 1, and are generally consistent with those reported in ref. ⁵ to within 1σ. However, >1σ differences are observed for a subset of parameters, namely: Δɸ for NRS1 is 3.7σ lower; Δɸ for NRS2 is 2.3σ lower; y_1,0 for NRS1 is 1.7σ lower; σ_syst for NRS1 is 1.5σ lower; and σ_syst for NRS2 is 3.3σ lower. These differences could arise from one or more of the analysis variations noted above. For example, the adoption of exoTiC-LD priors for the limb-darkening coefficients seems to have driven u₁ and u₂ away from the values found in ref. ⁵, which adopted substantially broader priors for those parameters. Nonetheless, the primary purpose of the white phase-curve fit was to refine the values of M_⋆, log₁₀R_⋆, M_p, ΔP_orb and i in particular, as these parameters were subsequently held fixed for the spectroscopic phase-curve fits (see next section). The updated values for the latter all agree to better than 0.8σ with those reported in ref. ⁵ and ref. ⁵⁶.

We also note that for the maximum a posteriori model—which in general we refer to as the ‘best-fit’ model—the root mean square (r.m.s.) of the residuals is found to be 126 ppm for the NRS1 phase curve and 153 ppm for the NRS2 phase curve, compared with Poisson uncertainties derived from the data reduction pipeline of 92 and 137 ppm, respectively. The scatter in the residuals bins down more slowly than expected for uncorrelated noise (Supplementary Fig. 2), which suggests that our instrument baseline model has not captured all of the systematics in the data. This should be inconsequential for the final atmospheric spectra that we obtain for WASP-121 b, as the white phase-curve fit is only used to refine our estimates for M_⋆, log₁₀R_⋆, M_p, ΔP_orb and i. Our posterior distributions for these parameters are fully consistent with the priors (Supplementary Fig. 1), which were adapted from the independent measurements cited above and in ref. ⁵. The sole exception is i, for which our refined value of ${{87.97}_{-0.15}^{+0.17}}^{\circ}$ is significantly lower than the prior value of 88.49 ± 0.16° taken from ref. ⁵². However, even if this ~0.5° difference between the prior and posterior for i is caused by uncorrected systematics in our white phase-curve fit, this would have a negligible effect on the final atmospheric spectra that we derive from the spectroscopic phase-curve fits described in the next section.

Spectroscopic phase-curve fitting

The spectroscopic phase curves were fitted in a similar fashion to the white phase curves. The main difference was that for each spectroscopic phase-curve fit, a subset of model parameters that are not expected to vary with wavelength were held fixed to the best-fit solution from the white phase-curve fit reported in Supplementary Table 1, namely: M_⋆, log₁₀R_⋆, M_p, ΔP_orb and i. All other parameters were allowed to vary separately in each spectroscopic phase-curve fit, including ΔT_mid and Δɸ. As for the white phase-curve fit, emcee was used to marginalize the posterior distribution for each spectroscopic phase curve, with 150 walkers, 500 burn-in steps and 5,000 sampling steps.

Supplementary Fig. 3 shows a selection of spectroscopic phase-curve fits spanning a range of fit qualities as measured by the χ² statistic. We find that the fit residuals are well-behaved and consistent with white noise. This is conveyed in the bottom panels of Supplementary Fig. 3, which provide a summary of the noise properties for the spectroscopic phase-curve fit residuals. The bottom left panel shows that the histograms of residuals taken across all wavelength channels are very close to normally distributed, with means close to zero and standard deviations within 2% of the pipeline Poisson noise. The bottom right panel shows the Allan deviations for the phase-curve residuals, with the r.m.s. of the residuals reducing at a rate of $\sqrt{N}$ where N is the number of data points per bin.

Generating the phase-resolved emission spectra

To extract the phase-resolved planetary emission spectra, we started by randomly selecting 5,000 posterior samples from each spectroscopic phase-curve fit. For each of these posterior samples, we generated the corresponding systematics baseline (β) and divided the associated data by this baseline to give a ‘corrected’ phase curve. We then binned the corrected data into 36 phase bins, each of 46 min duration. In defining the locations of the phase bins, we were careful to avoid the secondary eclipses and primary transit. Our chosen bins are illustrated in Extended Data Fig. 1, along with simplified representations of the relative fractions of the dayside and nightside hemispheres of the planet that were visible at each phase. The bin numbers are denoted n_ɸ, which increases in chronological order from n_ɸ = 1 for the first phase bin to n_ɸ = 36 for the last phase bin. In each phase bin, we took the median of the binned values generated from the 5,000 posterior samples as our measurement of F_p/F_⋆. To calculate the uncertainty, we took the quadrature sum $\sqrt{{\sigma }_{{\rm{r}}}^{2}+{\sigma }_{{\rm{s}}}^{2}}$, where σ_r denotes a random component and σ_s denotes a systematic component. For σ_r, we calculated the standard deviation of the model residuals in a given phase bin, then divided by ${\sqrt{N_j}}$, where N_j is the number of data points in the phase bin. We repeated this for all of the 5,000 posterior samples and took the median of these values as σ_r. To calculate σ_s, we first binned the phase-curve models (α) computed from each of the 5,000 posterior samples into the 36 phase bins. We then took the standard deviation of these binned α values in each phase bin as the σ_s value for that phase bin. This process was undertaken separately for each of the 349 wavelength channels to generate the full set of phase-resolved, wavelength-dependent F_p/F_⋆ measurements with uncertainties.

We repeated this process for independent sets of spectroscopic phase-curve fits performed using each of the instrument baselines (β₁, β₂ and β₃). We found that the resulting emission spectra obtained with β₁ and β₂ were consistent to within 1σ in all wavelength channels. The emission spectra obtained with β₁ and β₃ were consistent to within 1σ in all but one of the wavelength channels, with the values obtained for the latter channel still being consistent to within 1.2σ. β₁ was favoured by the BIC over β₂ and β₃ in 322 and 345 of the 349 wavelength channels, respectively. For the subset of channels in which either β₂ or β₃ was favoured by the BIC, the emission spectrum value obtained for each channel was always consistent to within 1σ of the corresponding value obtained with β₁. Given the overall BIC preference for β₁, along with the consistency of results obtained for each of the baselines tested, we adopt the emission spectra obtained using β₁ from this point onwards in our analysis.

Independent data reduction

As a check, we performed a second data reduction using the Eureka! pipeline⁵⁷. The steps we took were similar to those of ref. ⁵⁸, with a few modifications. For Stage 1, we set the up-the-ramp jump-detection threshold to 10σ, rather than the default 4σ, to avoid pixels being erroneously flagged as outliers. We corrected 1/f noise at the group-level within Eureka!’s Stage 1 by masking the trace and performing a column-by-column subtraction. For each column, we determined the offset to be subtracted by first excluding outliers >3 times the median of the unmasked pixels and fitting a zeroth-order polynomial to the remaining pixels in the column. We also applied a scale factor for the bias correction of each integration and both detectors using a smoothing filter calculated from the first group. Before extracting the spectrum from each integration, we corrected for the curvature of the G395H trace using the method described in ref. ⁵⁸ and performed a column-by-column background subtraction. For the background subtraction, we fitted a flat line to the pixels more than seven pixels away from the central row of the trace in that column. Before fitting to the background pixels, we performed two iterations of outlier masking, adopting outlier thresholds of 5σ along the time axis and >5 times the median along the spatial axis. The time-series spectra were then extracted from the background-subtracted integrations using optimal extraction⁵⁹ with an aperture of 7 pixels.

We generated phase curves from the Eureka! time-series spectra and fitted them using the same methodology described above for the FIREFLy data. Supplementary Fig. 4 compares the resulting nightside (n_ɸ = 18) and dayside (n_ɸ = 34) spectra to those obtained from the FIREFLy analysis. For the dayside and nightside, the pairs of spectra agree to within 1-σ for 344 and 341 of the 349 wavelength channels, respectively. Of the remaining channels, the agreement is better than 2σ for seven of the nightside channels and four of the dayside channels. This leaves a single channel of each spectrum for which the difference between the two analyses is >2σ. Namely, we obtain a 3.5σ difference for the 2.755 μm channel of the nightside spectrum, which is close to the short-wavelength edge of the NRS1 passband, and a similar 3.4σ difference for the 3.860 μm channel of the dayside spectrum. As these >2σ discrepancies are restricted to a single channel at each phase, they should not affect the interpretation of the spectra, which overall exhibit excellent agreement (Supplementary Fig. 4).

Atmospheric retrieval analyses

To characterize the dayside and nightside hemispheres of WASP-121 b, we performed retrieval analyses using five independent codes: ATMO^{60,61,62,63,64}, NEMESIS^65,66, CHIMERA⁶⁷, HyDRA^68,69,70 and PETRA⁷¹. The varying set-ups used for each of these codes are described below. Opacity sources considered by the retrieval codes^{72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108} are summarized in Supplementary Table 2. Further details of each retrieval code can be found in previous studies that have used them for exoplanet emission retrieval analyses, of which we cite a number of relevant examples here: ATMO was used in ref. ¹⁰⁹ and ref. ¹¹⁰ to analyse the dayside emission of WASP-121 b measured with HST and Spitzer; NEMESIS was used in ref. ⁵¹ to analyse the dayside and nightside spectra of WASP-43 b measured with JWST; CHIMERA was used in ref. ¹¹¹ to analyse the dayside and nightside spectra of the warm sub-Neptune GJ1214 b measured with JWST; HyDRA was used in ref. ¹¹² to analyse the dayside emission spectrum of the ultrahot giant planet WASP-18 b measured with JWST; and PETRA was used in ref. ⁷¹ to analyse the HST and Spitzer emission spectra of the hot giant planets HD 209458 b and WASP-43 b, as well as a synthetic JWST emission spectrum for the ultrahot giant planet KELT-9 b.

For this study of WASP-121 b, we have restricted our retrieval analyses to the n_ɸ = 34 and n_ɸ = 18 spectra, which are shown in Fig. 1. As can be seen in Extended Data Fig. 1, the n_ɸ = 34 spectrum corresponds to the phase bin immediately preceding the second secondary eclipse, while the n_ɸ = 18 spectrum corresponds to the phase bin immediately preceding the primary transit. In practice, we chose the n_ɸ = 34 and n_ɸ = 18 spectra to characterize the dayside and nightside hemispheres, respectively (Extended Data Fig. 1). However, it is important to note that retrieval results can potentially be biased if the measured emission at a given phase is modelled without allowing for inhomogeneous properties across the visible hemisphere^113,114. This is particularly relevant for the n_ɸ = 18 phase bin, where the nightside hemisphere dominates the fractional area coverage of the planetary disk, but the thin crescent of visible dayside hemisphere may still contribute meaningfully to the measured emission. We took two approaches to investigate how contamination by the dayside crescent affected the retrieval results for the n_ɸ = 18 spectrum. The first approach, taken for the NEMESIS retrieval, was to explicitly model each spectrum as a weighted sum of dayside and nightside contributions. The second approach, taken for the CHIMERA and HyDRA retrievals, was to include an additional free parameter in the form of the ‘dilution factor’ described in ref. ¹¹⁴.

The retrieval analyses presented in the following sections were performed on the 349-channel spectra (grey points in Fig. 1) and provided excellent fits to the data. Supplementary Fig. 5 shows histograms of the residuals for the best-fit models from each retrieval. In all cases, the residuals closely follow normal distributions. The reduced chi-squared values (${\chi }_{\upsilon }^{2}$) are also printed in each panel. For the ATMO equilibrium chemistry retrievals, we obtain ${\chi }_{\upsilon }^{2}=1.12$ for the dayside spectrum (n_ɸ = 34) and ${\chi }_{\upsilon }^{2}=1.09$ for the nightside spectrum (n_ɸ = 18). The ${\chi }_{\upsilon }^{2}$ values for all of the free chemistry retrievals (NEMESIS, PETRA, HyDRA and CHIMERA) fall between 1.04 and 1.10. Unless stated otherwise, reported results are median values with 1σ uncertainties defined by the 16th and 84th percentiles (68% credible interval).

ATMO retrieval

We used ATMO to perform independent retrievals for the n_ɸ = 18 and n_ɸ = 34 spectra. For each of the two retrievals, a single PT profile was assumed. Specifically, we adopted three-parameter¹¹⁵ and five-parameter^116,117,118 analytical forms to model the PT profiles for the n_ɸ = 18 and n_ɸ = 34 spectra, respectively. The five-parameter PT profile was adopted for the n_ɸ = 34 spectrum to allow more flexibility in capturing the shape of the dayside thermal inversion. We fixed the internal temperature to 500 K based on the results of ref. ¹¹⁹. By adopting a single PT profile for each phase, we did not explicitly allow for the possibility of temperature inhomogeneities across the visible hemisphere. However, as described below, additional tests performed with NEMESIS, CHIMERA and HyDRA demonstrate that this should not notably bias the final results.

For the chemistry, we allowed the carbon, oxygen and silicon abundances to vary as separate free parameters, while the relative abundances of all other heavy elements were held fixed to the solar values and varied jointly via an additional ‘metallicity’ parameter. ATMO adopts the solar abundances of ref. ¹²⁰, except for C, N, O, P, S, K and Fe, which are taken from ref. ³⁰. Given the PT profile and elemental abundances, the abundances of chemical species were then calculated assuming thermochemical equilibrium, which was done by minimizing the Gibbs free energy. The adopted priors for each model parameter are listed in Supplementary Table 3. The posterior distribution was sampled using nested sampling¹²¹.

The results of the ATMO retrieval are shown in Figs. 1–3. As already noted in the main text, spectral features of H₂O, CO, SiO and CH₄ allow the carbon, oxygen and silicon abundances to be constrained, as well as the atmospheric PT profile (Fig. 1). We confirm the dayside thermal inversion and infer a nightside PT profile that cools with decreasing pressure (Fig. 2), in agreement with previous studies of WASP-121 b^10,12. For the dayside (n_ɸ = 34), we obtain VMRs of ${\rm{C}}/{\rm{H}}={\mathrm{8,036}}_{-1,442}^{+1,771}\,{\mathrm{ppm}}$ and ${\rm{O}}/{\rm{H}}={\mathrm{8,741}}_{-1,516}^{+1,897}\,{\mathrm{ppm}}$, and derive a tight constraint for the atmospheric C/O ratio of ${0.92}_{-0.03}^{+0.02}$ (Fig. 3). We also obtain ${\rm{Si}}/{\rm{H}}={491}_{-138}^{+244}\,{\mathrm{ppm}}$ for the dayside, which, combined with the C/H and O/H measurements, allow us to derive ${\rm{C}}/{\rm{Si}}={15.78}_{-4.69}^{+7.22}$, ${\rm{O}}/{\rm{Si}}={17.07}_{-4.74}^{+7.49}$ and $({\rm{C}}+{\rm{O}})/{\rm{Si}}={32.84}_{-9.42}^{+14.73}$. For the nightside (n_ɸ = 18), we obtain ${\rm{C}}/{\rm{H}}={\mathrm{6,794}}_{-5,083}^{+10,399}\,{\mathrm{ppm}}$, which is consistent with the dayside constraint, but relatively weak owing to the lower signal-to-noise of the nightside emission. As described below for the free chemistry NEMESIS, PETRA, HyDRA and CHIMERA retrievals, we do not detect any O-bearing or Si-bearing species on the nightside, preventing us from obtaining useful constraints on the absolute O/H and Si/H abundances. However, we are able to constrain the relative nightside abundances to be C/O > 0.72 and C/Si > 0.74 at 95% probability, which are consistent with the dayside results.

The pressure-dependent VMRs for H₂O, CO, CO₂, SiO and CH₄ derived under the assumption of thermochemical equilibrium are shown for the nightside and dayside atmospheres in Fig. 2. At a pressure of 10⁻² bar on the nightside, we obtain base-10 logarithmic VMRs of ${-2.25}_{-0.44}^{+0.14}\,{\mathrm{dex}}$ for CH₄ and ${-2.73}_{-0.73}^{+0.29}\,{\mathrm{dex}}$ for CO. The equivalent logarithmic VMRs at a pressure of 10⁻² bar on the dayside are ${-12.38}_{-0.12}^{+0.10}\,{\mathrm{dex}}$ for CH₄ and ${-2.00}_{-0.10}^{+0.10}\,{\mathrm{dex}}$ for CO. For H₂O, CO₂ and SiO, thermal dissociation reduces the dayside abundances at pressures below approximately 10⁻² bar for H₂O and CO₂, and below approximately 10⁻³ bar for SiO. Note that photodissociation—which was not modelled—may further reduce molecular abundances at pressures below 10⁻⁵ bar (ref. ¹²²). Deeper in the dayside atmosphere at a pressure of 10⁻¹ bar, where neither thermal dissociation nor photodissociation reduce the abundances, the logarithmic VMRs are ${-3.57}_{-0.16}^{+0.18}\,{\mathrm{dex}}$ for H₂O, ${-6.13}_{-0.25}^{+0.24}\,{\mathrm{dex}}$ for CO₂ and ${-3.12}_{-0.14}^{+0.18}\,{\mathrm{dex}}$ for SiO. The nightside SiO abundance is weakly constrained, but nonetheless consistent with the inferred dayside abundance, while the nightside abundances of H₂O and CO₂ are found to be less than 1 ppm and 1 ppb, respectively.

Although abundances for many other species are predicted by the ATMO retrievals under the assumption of thermochemical equilibrium, we focus on H₂O, CO, CO₂, SiO and CH₄ because (with the exception of CO₂) they are the only ones that the free chemistry retrievals find direct evidence for in the data (see below). The detection of these molecules is also supported by the wavelength-dependent contribution functions for the ATMO retrievals, which are shown in Extended Data Fig. 3. The dayside contribution function traces out spectral features of H₂O, SiO and CO, while the nightside contribution function exhibits the distinctive spectral signature of CH₄. Supplementary Fig. 6 highlights the dayside SiO emission signal in particular, as it is relatively subtle compared with the dayside H₂O and CO signals. We discuss the nightside CH₄ signal in more detail below. From the contribution functions, we also note that the data seem to be probing a greater range of pressures on the nightside than the dayside. The absolute gradient inferred for the nightside temperature profile is much lower than that inferred for the dayside (Fig. 3), meaning that a relatively large range of pressures must be sampled to reproduce the spectral features observed in the nightside spectrum.

Extended Data Fig. 5 compares our retrieved PT profiles to those of various models that self-consistently treat the radiative transfer and chemistry of the atmosphere. Specifically, we show PT profiles for a range of dayside and nightside longitudes from published MITgcm¹¹ and Exo-FMS¹⁹ 3D GCM simulations, as well as 1D radiative-convective thermochemical equilibrium models computed with ATMO for internal temperatures of 100 K and 500 K, assuming the elemental abundances derived from the retrieval. Although the 3D models adopted solar elemental abundances, it is still interesting to make a preliminary comparison between the thermal profiles that they predict and those that we infer from the data. On the dayside, the retrieved PT profile is broadly consistent with the temperatures spanned by the self-consistent 1D and 3D models across the ~10⁻¹ to ~10⁻³ bar pressure range probed by the data (Fig. 2). On the nightside, the retrieved temperatures are similar to those predicted by the MITgcm simulation for pressures below ~10⁻² bar, but cooler at higher pressures. The Exo-FMS simulation predicts considerably warmer nightside temperatures at all pressures compared with our retrieved temperatures.

NEMESIS retrieval

We used NEMESIS to perform a joint retrieval for the n_ɸ = 18 and n_ɸ = 34 spectra. We considered the measured F_p at each phase to be the weighted sum of dayside and nightside components:

$${F}_{\mathrm{p}}=w{F}_{{\rm{day}}}+\left(1-w\right){F}_{{\rm{night}}},$$

(1)

where F_day and F_night denote the emission from the dayside and nightside hemispheres, respectively, and w is a geometric term giving the fraction of the visible hemisphere at orbital phase ɸ taken up by the dayside hemisphere¹²³. For the latter, we used:

$$w\left(\phi \right)=\left[1-\cos \left(2\uppi \phi -\Delta \phi \right)\right]/2$$

(2)

and fixed Δϕ = −3°, namely, the measured phase offset for the planetary brightness map rounded to the nearest degree (Supplementary Table 1 and Supplementary Fig. 1). We also allowed for the difference in the average emission angle from the dayside segment between the two phases. This could be particularly important for the dayside segment of the n_ɸ = 18 spectrum, which will be limb-darkened or limb-brightened in different spectral regions depending on the PT profile. We calculate the weighted-average emission angle of the dayside segment to be 76° at n_ɸ = 18 and 45° at phase n_ɸ = 34, and took these differences into account when computing the measured emission at each phase according to equation (1). We did not consider a similar correction for the nightside segment of phase n_ɸ = 34 to be necessary, because it makes a very small contribution to the overall spectrum measured at that phase.

We adopted a free chemistry approach for the NEMESIS retrieval, in which abundances of select chemical species were allowed to vary as free parameters without requiring thermochemical equilibrium to be satisfied. The chemical species considered were H₂O, H⁻, CO, CO₂, CH₄ and SiO. Abundances for these species were varied separately for the F_day and F_night model components. We allowed for the possibility of H₂O being thermally dissociated on the dayside by varying its abundance with P according to P^τ above some base pressure level P₀, which has been shown to provide a reasonable approximation for the effects of thermal dissociation¹¹. In addition to the H₂O abundance for P > P₀, we varied τ and P₀ as free parameters in the retrieval. All other abundances—including the nightside H₂O abundance—were assumed to be uniform in pressure. We adopted separate six-parameter PT profiles¹²⁴ for the dayside and nightside. Priors for each model parameter are listed in Supplementary Table 4. Posterior sampling was performed using the pymultinest¹²⁵ implementation of MultiNest^121,126 nested sampling.

As with the ATMO retrieval, the NEMESIS retrieval recovers a thermal inversion on the dayside and a temperature profile that cools with decreasing pressure on the nightside (Supplementary Fig. 7). Tight abundance constraints indicate the presence of H₂O, CO and SiO on the dayside, and CH₄ on the nightside (Supplementary Fig. 7). We infer logarithmic VMRs of ${-3.31}_{-0.26}^{+0.23}\,{\mathrm{dex}}$ for the deep-atmosphere H₂O, ${-1.29}_{-0.19}^{+0.13}\,{\mathrm{dex}}$ for CO and ${-2.66}_{-0.21}^{+0.16}\,{\mathrm{dex}}$ for SiO on the dayside and ${-1.75}_{-0.66}^{+0.44}\,{\mathrm{dex}}$ for CH₄ on the nightside. We also find that the inferred dayside and nightside spectra extrapolate well to the other orbital phases (Extended Data Fig. 2).

To evaluate the significance of the CH₄ detection, we performed two nightside-only NEMESIS retrievals, with CH₄ included and with CH₄ removed. In these additional retrievals, we also included C₂H₂, HCN and NH₃, to ensure that the signal we observe really is due to CH₄. By comparing the Bayesian evidence for this ‘null hypothesis’ retrieval without CH₄ (${{\mathcal{Z}}}_{0}$) to the Bayesian evidence of the model that included all of the molecules (${{\mathcal{Z}}}_{1}$), we obtain a log Bayes factor ($\mathrm{ln}{{\mathcal{B}}}_{10}=\mathrm{ln}{{\mathcal{Z}}}_{1}/{{\mathcal{Z}}}_{0}$) of $\mathrm{ln}{{\mathcal{B}}}_{10}=10.4$. Under numerous empirically calibrated evidence scales, this corresponds to the highest band of evidence in favour of a CH₄ detection^127,128,129. Following refs. ^130,131,132, we can also translate our $\mathrm{ln}{{\mathcal{B}}}_{10}$ value into an equivalent frequentist detection significance, as is common in the exoplanet literature. Using this approach, we obtain a detection significance of 4.9σ for CH₄ on the nightside of WASP-121 b. Similar dayside-only retrievals were performed to obtain decisive detections of H₂O ($\mathrm{ln}{{\mathcal{B}}}_{10}=88.3$, 13.5σ), CO ($\mathrm{ln}{{\mathcal{B}}}_{10}=56.0$, 10.8σ) and SiO ($\mathrm{ln}{{\mathcal{B}}}_{10}=14.5$, 5.7σ). Plots of the best-fit NEMESIS retrievals with and without each molecule are shown in Supplementary Figs. 8–11.

We additionally used NEMESIS to fit for a nightside cloud layer using the cloud parameterization of ref. ¹³³, as was done previously in ref. ¹³⁴. We fitted for the top pressure of an opaque, grey cloud (P_top), and a scattering power-law index (γ) and opacity scaling (σ₀) for a ‘haze’ above the cloud. From this retrieval, we are only able to provide a weak constraint for the cloud-top pressure being >10⁻³ bar, while the remaining cloud parameters are effectively unconstrained.

CHIMERA retrieval

We used CHIMERA to perform free chemistry retrievals for the n_ɸ = 18 spectrum, focusing on characterization of the nightside atmosphere only. The chemical species considered were H₂O, CO, CO₂, CH₄, C₂H₂, HCN and NH₃. We assumed that all VMRs were uniform with pressure, as the temperatures on the nightside are too low for thermal dissociation to be important. We performed separate retrievals adopting widely used three-parameter¹¹⁵, five-parameter¹¹⁶ and six-parameter¹²⁴ functions for the PT profile. The three-parameter profile was the same as that used for the ATMO nightside analysis, the five-parameter profile was the same as that used for the ATMO dayside analysis and the six-parameter profile was the same as that used for the NEMESIS analyses. We also performed retrievals with and without the dilution factor described in ref. ¹¹⁴ to account for contribution to the observed spectrum from the visible dayside crescent. The priors used for each model parameter are listed in Supplementary Table 5. Posterior samples were drawn using pymultinest¹²⁵.

The retrievals performed for each PT profile support the presence of CH₄, with $\mathrm{ln}{{\mathcal{B}}}_{10}=5.60$ (3.8σ) for the three-parameter profile, $\mathrm{ln}{{\mathcal{B}}}_{10}=9.32$ (4.7σ) for the five-parameter profile and $\mathrm{ln}{{\mathcal{B}}}_{10}=11.0$ (5.1σ) for the six-parameter profile. The inferred PT profiles are plotted in Supplementary Fig. 7 and all cool with decreasing pressure across the 1–10⁻³ bar pressure range. The six-parameter profile is essentially identical to the NEMESIS six-parameter profile and is also in good agreement with the HyDRA six-parameter profile (see below). The three-parameter profile is fully consistent with the ATMO three-parameter profile. Interestingly, the five-parameter profile exhibits an inversion for pressures below approximately 10⁻³ bar, producing a narrow CH₄ emission spike at 3.35 μm while the broader CH₄ features remain in absorption, which marginally improves the fit to the data (Supplementary Fig. 12). However, on the basis of the evaluated Bayesian evidence values, we do not find a clear preference for any one of the profiles in particular. Specifically, we obtain $\mathrm{ln}{{\mathcal{Z}}}_{1}=-192.72\pm 0.07$ for the three-parameter profile, $\mathrm{ln}{{\mathcal{Z}}}_{1}=-190.39\pm 0.08$ for the five-parameter profile and $\mathrm{ln}{{\mathcal{Z}}}_{1}=-191.56\pm 0.06$ for the six-parameter profile. A comparison of the CHIMERA best-fit models with and without CH₄ for the five-parameter profile is shown in Supplementary Fig. 12.

As shown in Supplementary Fig. 7, the corresponding logarithmic VMRs for CH₄ are found to be ${-2.43}_{-1.67}^{+0.91}\,{\mathrm{dex}}$ (three-parameter), ${-3.73}_{-0.57}^{+0.49}\,{\mathrm{dex}}$ (five-parameter) and ${-1.73}_{-0.78}^{+0.48}\,{\mathrm{dex}}$ (six-parameter). Besides CH₄, no evidence is obtained for any of the other chemical species that were considered. We also find that the inclusion of a dilution factor does not significantly affect the results and is not justified by the Bayesian evidence.

HyDRA retrieval

We used HyDRA to perform separate free chemistry retrievals for the n_ɸ = 18 and n_ɸ = 34 spectra. The chemical species considered for the dayside (n_ɸ = 34) were H₂O, CO, CO₂, SiO, SiO₂, HCN, H⁻, OH, FeH, TiO and VO. For the nightside (n_ɸ = 18), we considered H₂O, CO, CO₂, SiO, SiO₂, HCN, H⁻ and CH₄. Thermal dissociation was modelled for the dayside abundances of H₂O, TiO, VO and H⁻ using the methods of ref. ¹¹. The deep abundances of each of these species were free parameters in the retrieval, while their pressure- and temperature-dependent abundance profiles were calculated using the coefficients in table 1 of ref. ¹¹ as inputs to equations (1) and (2) of the same work. All other VMRs were assumed to be uniform in pressure. The six-parameter function¹²⁴ was used to model the PT profiles at both phases. We tested the inclusion of clouds in the retrievals, which were parameterized by a modal particle size (a_c), the cloud base pressure (P_c), a pressure exponent (ε), the cloud particle abundance (f_c) and the cloud covering fraction (ψ_c). The cloud particle abundance was assumed to be zero at pressures greater than P_c, and to decrease at pressures less than P_c, such that at pressure P the abundance was f_c(P/P_c)^ε. We assumed an MgSiO₃ cloud composition, and used the cloud absorption cross-sections from ref. ¹³⁵. The priors adopted for each model parameter are listed in Supplementary Table 6. Posterior sampling was performed using pymultinest¹²⁵.

The HyDRA retrieval again recovers a thermal inversion for the dayside and a cooling PT profile for the nightside. On the dayside, decisive detections are made for CO with $\mathrm{ln}{{\mathcal{B}}}_{10}=78.7$ (12.8σ), SiO with $\mathrm{ln}{{\mathcal{B}}}_{10}=17.5$ (6.2σ) and H₂O with $\mathrm{ln}{{\mathcal{B}}}_{10}=13.4$ (5.5σ). The inferred logarithmic VMRs for CO and SiO on the dayside are $-{2.17}_{-0.37}^{+0.35}\,{\mathrm{dex}}$ and $-{3.17}_{-0.34}^{+0.33}\,{\mathrm{dex}}$, respectively, while the logarithmic VMR for the deep atmosphere H₂O is found to be $-{4.24}_{-0.32}^{+0.35}\,{\mathrm{dex}}$. We caution that the inferred SiO abundance may be biased, as it was assumed to be uniform in pressure, but in reality should be thermally dissociated (Fig. 2). On the nightside, the detection significance for CH₄ is found to be $\mathrm{ln}{{\mathcal{B}}}_{10}=3.3$ (3.1σ), with an inferred logarithmic VMR of ${-2.05}_{-0.81}^{+0.55}\,{\mathrm{dex}}$. However, we also ran a separate retrieval analysis for the nightside spectrum with CH₄ as the only opacity source and compared the Bayesian evidence to that obtained for a simple blackbody fit. In this restricted scenario, HyDRA decisively favours a model with CH₄ absorption over a blackbody with $\mathrm{ln}{B}_{10}=31.5$ (8.1σ). Plots of the best-fit HyDRA retrievals with and without each detected molecule are shown in Supplementary Figs. 13–16. We additionally find that a dilution parameter is not statistically justified in either the dayside or nightside retrievals, nor are clouds.

PETRA retrieval

We used PETRA to perform free chemistry retrievals for the n_ɸ = 34 spectrum, focusing on characterization of the dayside atmosphere only. The chemical species we considered were H₂O, CO and SiO. Only the H₂O abundance was allowed to be affected by thermal dissociation, following a similar approach to that used in the NEMESIS and HyDRA free chemistry retrievals. The PT profile was modelled with the same five-parameter function used in the ATMO, CHIMERA and HyDRA retrievals. The priors used for each model parameter are listed in Supplementary Table 7. Posterior sampling was performed using differential-evolution Markov chain Monte Carlo^136,137.

As for all other dayside retrievals, the inferred PT profile exhibits a strong thermal inversion (Supplementary Fig. 7). However, the thermal inversion inferred by PETRA starts at a pressure of ~10⁻³ bar, which is higher in the atmosphere than the thermal inversions inferred by the other retrievals. For the chemical abundances, we obtain logarithmic VMRs of ${-2.03}_{-0.36}^{+0.32}\,{\mathrm{dex}}$ for CO, ${-3.15}_{-0.28}^{+0.29}\,{\mathrm{dex}}$ for SiO and ${-3.92}_{-0.27}^{+0.28}\,{\mathrm{dex}}$ for H₂O. We note that the lower edge of the 68% credible interval (16th percentile) for the H₂O abundance decreases with increasing pressure for pressures above approximately 10⁻² bar. This corresponds to solutions where the thermal dissociation exponent is negative (τ < 0) and indicates that the data are not constraining the H₂O abundance for pressures >10⁻² bar in the PETRA retrieval.

Discussion of the retrieval analyses

The four retrievals that considered the n_ɸ = 34 spectrum (ATMO, NEMESIS, HyDRA and PETRA) recover dayside PT profiles with thermal inversions and fit the data with emission features of H₂O, CO and SiO. The exact shapes of the inferred PT profiles differ somewhat and appear to be correlated with differences in the inferred pressure-dependent abundances, particularly that of H₂O. As shown in Supplementary Fig. 7, NEMESIS and PETRA favour the dissociation of H₂O at pressures <10⁻³ bar, resulting in thermal inversions at lower pressures than those obtained by ATMO and HyDRA, for which H₂O dissociation is found to occur at pressures of ~10⁻² bar. These results emphasize how challenging it can be to model pressure-variable abundances with free chemistry retrievals in a purely parametric manner, even when constrained by JWST-quality data. Despite the variation in absolute abundances and thermal inversion pressures, the inferred C/O, C/Si and O/Si ratios are found to be highly consistent across all the retrieval analyses (Supplementary Fig. 17).

Each of the retrieval analyses that were performed for the n_ɸ = 18 spectrum (ATMO, NEMESIS, HyDRA and CHIMERA) return nightside PT profiles that cool with decreasing pressure and find evidence for CH₄ absorption features. The PT profiles and absolute CH₄ abundances are in good overall agreement across the retrievals (Supplementary Fig. 7). The one mild exception is the CHIMERA retrieval that adopted the five-parameter PT profile, which finds a low-pressure inversion that is not seen in the other retrieval analyses. As noted above, because this inversion occurs at low pressures it only has a minor effect on the appearance of the emission spectrum. Otherwise, the CHIMERA five-parameter PT profile broadly resembles those of the other analyses, but shifted to higher pressures. This reflects the relatively low CH₄ abundance inferred by CHIMERA for the five-parameter PT profile (Supplementary Fig. 7), which has the effect of moving the photosphere to a higher pressure than for the other retrieval models.

In assessing the detection significances for each molecule, we report the corresponding range of values obtained from the various free chemistry retrievals: 5.5–13.5σ for dayside H₂O; 10.8–12.8σ for dayside CO; 5.7–6.2σ for dayside SiO; and 3.1–5.1σ for nightside CH₄. Within these ranges, we favour those detection significances obtained by the retrieval codes that used the latest line list for the corresponding molecule (Supplementary Table 2). For H₂O, HyDRA uses the HITEMP line list of ref. ⁷⁴, which is based on the earlier ab initio line list of ref. ⁷², while both NEMESIS and CHIMERA use the more recent ExoMOL line list of ref. ⁷³. For CO, NEMESIS uses the ExoMOL line list of ref. ⁷⁹, while HyDRA and CHIMERA use the earlier HITEMP line list of ref. ⁷⁴. For SiO, HyDRA uses the latest ExoMOL line list of ref. ⁸², which superseded the earlier line list of ref. ⁸¹ used by NEMESIS. For CH₄, CHIMERA uses the recent HITEMP line list of ref. ⁷⁸, whereas NEMESIS and HyDRA use the earlier ExoMOL line lists of refs. ^75,76,77. On the basis of these line lists, our favoured detection significances are therefore 13.5σ for dayside H₂O (NEMESIS), 10.8σ for dayside CO (NEMESIS) and 6.2σ for dayside SiO (HyDRA). For nightside CH₄, we favour the detection significance of 4.7σ from the CHIMERA retrieval performed with the five-parameter PT profile, which gave a higher Bayesian evidence ($\mathrm{ln}{{\mathcal{Z}}}_{1}=-190.39\pm 0.08$) than the three-parameter ($\mathrm{ln}{{\mathcal{Z}}}_{1}=-192.72\pm 0.07$) and six-parameter ($\mathrm{ln}{{\mathcal{Z}}}_{1}=-191.56\pm 0.06$) PT profiles.

We note that the lowest detection significance of 3.1σ obtained for nightside CH₄ by the HyDRA retrieval may be overly conservative. Higher detection significances of 5.1σ and 4.9σ, respectively, were obtained by the CHIMERA and NEMESIS retrievals for the same six-parameter PT profile adopted by HyDRA. We suspect that the difference in detection significances may be due to HyDRA using the early ExoMOL line lists of refs. ^75,77 for CH₄, whereas NEMESIS uses the expanded ExoMOL line list of ref. ⁷⁶, and CHIMERA uses the more recent HITEMP line list of ref. ⁷⁸, which (unlike the ab initio ExoMOL line lists) is largely based on experimental data. Differences in the adopted H₂O line list could also potentially affect the inferred CH₄ detection significance, given that H₂O and CH₄ have broad features at similar wavelengths within the G395H passband. As noted above, NEMESIS and CHIMERA both use the H₂O line list of ref. ⁷³, which is more recent than the line list of ref. ⁷⁴ used by HyDRA. Indeed, the dayside H₂O detection significance obtained by NEMESIS (13.5σ) is considerably higher than that obtained by HyDRA (5.5σ), as would be expected if the H₂O line list used by NEMESIS is better suited for the conditions of the WASP-121 b atmosphere than the older line list used by HyDRA. Furthermore, the fact that NEMESIS obtains lower detection significances than HyDRA for the other molecules (that is, CO and SiO on the dayside) demonstrates that NEMESIS does not systematically overestimate detection significances relative to HyDRA.

In addition to the formal detection significances, there are several other reasons that we favour nightside models including CH₄. First, when CH₄ is removed from the free chemistry retrievals (NEMESIS, HyDRA and CHIMERA), a thermal inversion at pressures >10⁻³ bar is required to reproduce the observed nightside spectrum with emission features of H₂O and CO. A thermal inversion this deep in the non-irradiated nightside atmosphere is not expected, on the basis of state-of-the-art GCM simulations (Extended Data Fig. 5). It would also be in tension with HST measurements that have previously revealed a decreasing PT profile on the nightside of WASP-121 b¹². Second, the inference of an atmospheric ${\rm{C}}/{\rm{O}}={0.92}_{-0.03}^{+0.02}$ from the dayside spectrum (ATMO) is a robust result (Fig. 3). With an atmospheric C/O ratio this high, equilibrium chemistry predicts a logarithmic VMR above −4 dex for CH₄ at the temperatures of the nightside (Extended Data Fig. 4), which is broadly consistent with the abundances inferred by the free chemistry retrievals (Supplementary Fig. 7). Although horizontal quenching of CH₄-poor gas from the dayside could in principle dominate the nightside composition, this nonetheless demonstrates that the CH₄ abundances we have inferred are physically plausible. We also reiterate that we did allow for dayside contamination of the n_ɸ = 18 spectrum using two approaches: including a dilution factor as a free parameter (HyDRA and CHIMERA) and by explicitly modelling the dayside contribution (NEMESIS). The evidence for CH₄ and a decreasing PT profile on the nightside persists with both of these approaches.

We do not make a statistically significant detection of CO on the nightside, despite equilibrium chemistry predicting a CO abundance that is comparable to the CH₄ abundance (Fig. 2 and Extended Data Fig. 4). This is because the abundance-weighted CO opacity is only higher than that of CH₄ within narrow line cores across the ~4–4.5 μm wavelength range. At the spectral resolution of the nightside data, the corresponding opacity enhancement provided by CO over CH₄ is quite modest (Extended Data Fig. 3) and we find that the inclusion or exclusion of CO has only a minor effect on the spectrum relative to the measurement uncertainties (Fig. 1). In particular, when we simply switch off the CO opacity in the best-fit ATMO model shown in Fig. 1, the χ² increases from 370.6 to 372.4, corresponding to a negligible ${\chi }_{\upsilon }^{2}$ increase from 1.087 to 1.092. Consequently, the nightside CO abundance is only weakly constrained by the free chemistry retrievals, while the CO abundance displayed in Fig. 2 primarily reflects the chemical equilibrium constraint imposed for the ATMO retrieval.

We also do not find evidence for clouds on the nightside of WASP-121 b, based on the Bayes factor computed from the HyDRA retrievals performed with and without clouds, as well as the lack of cloud constraints obtained by NEMESIS. This is despite the inferred nightside PT profile crossing numerous condensation curves of refractory species (Fig. 2), including silicates, which are expected to be the dominant source of cloud opacity on the hot nightsides of planets like WASP-121 b¹⁴. However, as noted in the main text, our non-detection of clouds does not necessarily imply the absence of nightside clouds. In particular, the inferred PT profile for the nightside hemisphere crosses the MgSiO₃, Mg₂SiO₄ and SiO₂ condensation curves at pressures of ~1 bar, whereas the data are only sensitive to pressures below ~10⁻¹ bar (Fig. 2 and Extended Data Fig. 6). It therefore remains possible that a thick cloud layer composed of silicates and/or other refractory species, such as Fe, is present on the nightside of WASP-121 b at pressures >10⁻¹ bar.

Stellar abundances and relative planetary abundances

We derived photospheric and fundamental stellar parameters for the WASP-121 host star using the algorithm described in ref. ¹³⁸. For our isochrone fitting, we included multiwavelength photometry from the UV to the near-infrared: Tycho-2 BT and VT¹³⁹, Gaia Data Release 2 (DR2)^{140,141,142,143}, the G, J, H and Ks bands from the Two Micron All Sky Survey All-Sky Point Source Catalog¹⁴⁴ and the W1 and W2 bands from the Wide-field Infrared Survey Explorer AllWISE mid-infrared data^145,146. We also included the Gaia DR3¹⁴³ parallax-based distance to the WASP-121 system¹⁴⁷. We included the extinction A_V inference based on 3D maps of extinction in the solar neighbourhood from the Structuring by Inversion the Local Interstellar Medium programme^148,149,150.

For the spectroscopic-based inferences, we used the equivalent widths of Fe i and Fe ii atomic absorption lines. The equivalent widths were derived from a spectrum of the WASP-121 host star measured with the ESPRESSO instrument on the Very Large Telescope. The spectrum is available from the European Southern Observatory archive under programme ID 106.21QM.001. It was taken on 2021 September 6 and has a signal-to-noise ratio of ~150 at 5,000 Å. To measure the equivalent widths, we used the splot task of IRAF (IRAF is distributed by the National Optical Astronomy Observatory, which is operated by the Association of the Universities for Research in Astronomy, Inc. under a cooperative agreement with the National Science Foundation). The atomic data and measured equivalent widths for each line are reported in Supplementary Table 8. We adapted the atomic data from ref. ¹⁵¹ and assumed photospheric solar abundances from ref. ¹⁵².

We used the isochrones package¹⁵³ to fit the MESA Isochrones and Stellar Tracks^154,155,156 library to our photospheric stellar parameters, as well as our input multiwavelength photometry, parallax and extinction data using pymultinest¹²⁵. Further details of our approach are provided in ref. ¹³⁸. The adopted stellar parameters are given in Supplementary Table 9.

We inferred elemental abundances of C i, O i, Na i, Mg i, Al i, Si i, S i, K i, Ca i, Sc ii, Ti ii, V i, Cr i, Cr ii, Mn i, Fe i, Fe ii, Co i, Ni i, Cu i, Zn ii, Y ii, Zr ii, Ba ii, Ce ii, Nd ii and Sm ii using the equivalent-width method, including isotopic/hyperfine splitting details where needed. We measured the equivalent widths by fitting Gaussian profiles with the splot task in IRAF. We also used the deblend task to disentangle absorption lines from adjacent spectral features whenever necessary. We assumed the solar abundances of ref. ¹⁵² and local thermodynamic equilibrium (LTE) and used the 1D plane-parallel solar-composition MARCS model atmospheres¹⁵⁷ and the 2019 version of the LTE radiative transfer code MOOG¹⁵⁸ to infer elemental abundances based on our equivalent widths. Our abundance inferences and uncertainties are reported in Supplementary Table 10. Where possible, we updated our elemental abundances for departures from LTE by linearly interpolating published grids of non-LTE corrections. Specifically, non-LTE corrections taken from the tables of refs. ^{159,160,161,162} were applied to the carbon, oxygen, aluminium, calcium, silicon, potassium and iron abundances.

We highlight the stellar abundances derived for C i, O i and Si i in particular, with reference to the abundances in the standard A(X) format. Here A(X) = logN_X/N_H + 12, where N_X is the number density of species X and N_H is the number density of hydrogen. The non-LTE (LTE) abundances we obtain are: 8.527 ± 0.050 (8.545 ± 0.050) for C i; 8.855 ± 0.024 (9.238 ± 0.037) for O i; and 7.696 ± 0.097 (7.711 ± 0.097) for Si i. The non-LTE and LTE abundances are consistent to well within the 1σ uncertainties for both C i and Si i. However, the non-LTE abundance for O i is significantly lower than the LTE abundance. This results in an upward revision of the stellar C/O ratio from 0.203 ± 0.050 (LTE) to 0.470 ± 0.061 (non-LTE). As an aside, we note that our LTE C/O ratio is consistent with the value of 0.23 ± 0.05 reported in ref. ¹⁶³.

To compute elemental ratios for the planetary atmosphere relative to those of the host star, we used the ATMO posterior samples shown in Fig. 3 for the planet abundances and randomly drew stellar abundances from normal distributions that had means and standard deviations set to the non-LTE values listed in Supplementary Table 10. We find that the C/H, O/H, Si/H and C/O ratios of the planet normalized to the equivalent ratios of the host star are greater than the following values at >99.99% probability: (C/H)/(C/H)_⋆ > 11.92, (O/H)/(O/H)_⋆ > 7.27, (Si/H)/(Si/H)_⋆ > 2.49 and (C/O)/(C/O)_⋆ > 1.63. The corresponding medians with uncertainties reflecting the 68% credible intervals (16th–84th percentiles) are: ${({\rm{C}}/{\rm{H}})/({\rm{C}}/{\rm{H}})}_{\star }={23.96}_{-4.75}^{+6.13}$, ${({\rm{O}}/{\rm{H}})/({\rm{O}}/{\rm{H}})}_{\star }={12.19}_{-2.24}^{+2.78}$, ${({\rm{Si}}/{\rm{H}})/({\rm{Si}}/{\rm{H}})}_{\star }={9.89}_{-2.89}^{+5.99}$ and ${({\rm{C}}/{\rm{O}})/({\rm{C}}/{\rm{O}})}_{\star }$$={1.96}_{-0.11}^{+0.11}$. We also find that the volatile/refractory ratios of the planet normalized to the host star are less than the following values at >99.99% probability: (C/Si)/(C/Si)_⋆ < 8.07, (O/Si)/(O/Si)_⋆ < 4.02 and ((C + O)/Si)/((C + O)/Si)_⋆ < 5.24. The corresponding 68% credible intervals are ${({\rm{C}}/{\rm{Si}})/({\rm{C}}/{\rm{Si}})}_{\star }={2.30}_{-0.61}^{+1.06}$, ${({\rm{O}}/{\rm{Si}})/({\rm{O}}/{\rm{Si}})}_{\star }={1.19}_{-0.32}^{+0.51}$ and ${(({\rm{C}}+{\rm{O}})/{\rm{Si}})/(({\rm{C}}+{\rm{O}})/{\rm{Si}})}_{\star }={1.55}_{-0.41}^{+0.71}$.

Accretion of rocky material

We can use the Si/H enrichment of the planetary atmosphere relative to the host star to estimate the quantity of rocky material that WASP-121 b accreted following its formation. If we assume that the planet’s primordial atmospheric envelope has Si/H equal to that of the host star, and is then enriched by silicon through the subsequent accretion of rocky material, then:

$${n}_{{\mathrm{Si}},{\mathrm{env}}}/{n}_{{\mathrm{H}},{\mathrm{env}}}={n}_{{\mathrm{Si}},\star }/{n}_{{\mathrm{H}},\star }+{n}_{{\mathrm{Si}},\upvarepsilon }/{n}_{{\mathrm{H}},{\mathrm{env}}},$$

where n_Si,env is the total number of silicon atoms in the atmospheric envelope, n_H,env is the total number of hydrogen atoms in the atmospheric envelope, n_Si,⋆/n_H,⋆ is the ratio of silicon atoms to hydrogen atoms in the host star and n_Si,ε is the number of accreted silicon atoms. Note that the above formula implicitly assumes that the accreted silicon is mixed throughout the entire atmospheric envelope.

The number of hydrogen atoms in the atmospheric envelope is given by:

$${n}_{{\mathrm{H}},{\mathrm{env}}}=\frac{{N}_{\mathrm{A}}\;{f}_{\mathrm{H}}({M}_{\mathrm{p}}-{M}_{{\mathrm{core}}})}{{\mu }_{\mathrm{H}}}$$

where N_A is Avogadro’s constant, f_H is the hydrogen mass fraction for the envelope and μ_H is the molar mass of hydrogen. We have also assumed that M_p = M_core + M_env, where M_core is the mass of the core and M_env is the mass of the atmospheric envelope.

Combining the above, we derive the following expression for the number of accreted silicon atoms:

$${n}_{{\mathrm{Si}},\upvarepsilon }=\frac{{N}_{\mathrm{A}}\;{f}_{\mathrm{H}}({M}_{\mathrm{p}}-{M}_{{\mathrm{core}}})}{{\mu }_{\mathrm{H}}}\left[{\left(\frac{{n}_{{\mathrm{Si}}}}{{n}_{\mathrm{H}}}\right)}_{{\mathrm{env}}}-{\left(\frac{{n}_{{\mathrm{Si}}}}{{n}_{\mathrm{H}}}\right)}_{\star }\right].$$

To calculate n_Si,ε, we take ${\left({n}_{{\mathrm{Si}}}/{n}_{\mathrm{H}}\right)}_{{\mathrm{env}}}={491}_{-138}^{+244}\,{\mathrm{ppm}}$, ${\left({n}_{{\mathrm{Si}}}/{n}_{\mathrm{H}}\right)}_{\star }=50\,{\mathrm{ppm}}$ and ${M}_{\mathrm{p}}=1.2{M}_{\mathrm{J}}$ (where M_J is the mass of Jupiter) from our measurements reported above, as well as the standard values for N_A (6.022 × 10²³ atoms mol⁻¹) and μ_H (10⁻³ kg mol⁻¹). We also assume M_core = 15 M_⊕, based on the estimate of Jupiter’s core mass¹⁶⁴, and f_H = 75%, based on the approximate hydrogen mass fraction of Jupiter’s atmosphere¹⁶⁵. With these input values, we obtain n_Si,ε = 1.63 × 10⁵⁰ atoms of silicon added to the envelope of WASP-121 b by the accretion of rocky material.

If we further assume that the rocky material has an Earth-like composition, we can translate this number of accreted silicon atoms to an equivalent number of Earth masses of rocky material. The estimated silicon mass fraction for the bulk Earth is 16% (ref. ¹⁶⁶). Given that Earth has a mass of 5.9723 × 10²⁴ kg, this equates to 9.5557 × 10²³ kg of silicon in the Earth. Adopting μ_Si = 0.028 kg mol⁻¹ for the molar mass of silicon, the total number of silicon atoms in the Earth then works out to be approximately n_Si,⊕ = 2.06 × 10⁴⁹. Therefore, we estimate that ${n}_{{\mathrm{Si}},\upvarepsilon }/{n}_{{\mathrm{Si}},\bigoplus }={21.2}_{-6.6}^{+11.8}$, implying that the accretion of approximately 15–33 M_⊕ of rocky material with an Earth-like composition can explain the silicon enrichment we have measured for the atmosphere of WASP-121 b.

As noted above, this calculation assumes that the planetary envelope started out with a silicon enrichment equal to the host star. However, the gas accreted by the planet at formation may have been depleted in silicon, as silicon would have condensed out of the gas phase throughout the protoplanetary disk^2,24. If we repeat the above calculation assuming that all of the silicon observed in the atmosphere was acquired through the accretion of rocky material by setting n_Si,env/n_H,env = n_Si,ε/n_H,env, then we estimate the quantity of accreted rocky material to be ${n}_{{\mathrm{Si}},\upvarepsilon }/{n}_{{\mathrm{Si}},\bigoplus }={23.7}_{-6.6}^{+11.8}$, or approximately 17–36 M_⊕.

Other wavelength-dependent quantities of interest

Supplementary Fig. 18 shows the atmospheric transmission spectrum of the day–night terminator region given by the inferred values for R_p as a function of wavelength. The data exhibit a downward slope between 2.75 μm and 4 μm, then increase sharply and peak at around 4.5–4.8 μm, before dropping off again towards longer wavelengths. These wavelength-dependent variations in the transmission spectrum indicate that the data are precise enough to identify key absorbing species at the day–night terminator. Predictions from models published in ref. ¹ and ref. ¹⁶⁷ are also shown in Supplementary Fig. 18, both of which were fitted to optical and near-infrared data obtained with HST. Neither model provides a good fit to the longer-wavelength JWST data. Further analysis of the transmission spectrum will be presented in an accompanying paper¹⁶⁸.

Wavelength-dependent values for u₁ and u₂, ΔT_mid and Δɸ are also shown in Supplementary Figs. 19–21. We find that the posterior values for u₁ and u₂ closely match the exoTiC-LD priors. For ΔT_mid and Δɸ, we do not see any clear evidence for variations as a function of wavelength that can easily be explained. A flat line fitted to the ΔT_mid values gives ${\chi }_{\upsilon }^{2}=1.00$, as would be expected if there are no wavelength-dependent variations. However, a higher ${\chi }_{\upsilon }^{2}=1.74$ is obtained when a flat line is fitted to the Δɸ values, indicating possible evidence for wavelength-dependent variations that should be investigated further in future studies.