Main

We observed one transit of the exoplanet WASP-107b (ref. 12) with the G395H spectral grating of the James Webb Telescope Near-Infrared Spectrograph (JWST-NIRSpec) as part of Guaranteed Time Observations (GTO) programme 1224 (principal investigator S. Birkmann). In the first years of operation of the James Webb Space Telescope (JWST), this mode demonstrates reliable detections of H2O, CO, CO2 and SO2 for several giant exoplanets7,13,14,15. See the Methods for further details on the observations and data analysis. The wavelength-integrated JWST-NIRSpec time-series photometry of WASP-107b is shown in Fig. 1, and the transmission spectrum is shown in Fig. 2. The spectrum shows several molecular absorption signals that are readily identifiable. We detect a large CO2 absorption feature at 4.3 μm (ref. 15) (Fig. 2) and a slope-like feature at the shortest wavelengths from H2O (ref. 14).

Fig. 1: The light curve of WASP-107b observed by JWST-NIRSpec G395H.
figure 1

a, The normalized wavelength-integrated white-light curves for the two detectors are shown, with the NRS1 (2.70–3.71 μm) and NRS2 (3.83–5.16 μm) detectors offset for clarity. A best-fit limb-darkened transit model is overplotted (blue). See Extended Data Fig. 2 for further details. b, The model residuals that achieve near-photon limited precisions.

Source Data

Full size image
Fig. 2: WASP-107b transmission spectral measurements.
figure 2

JWST-NIRSpec transmission spectrum and the 1σ uncertainties. The best-fit ATMO32 model is also plotted, and the individual contributions for each molecular species are shown.

Source Data

Full size image

We ran a series of model retrievals to provide detailed constraints on the atmospheric composition and temperature (Methods). The H2O and CO2 features are detected at high confidence (Extended Data Table 2). The spectrum shows CO features between 4.5 μm and 5 μm (4σ), and an SO2 absorption feature at 4.0 μm (5.5σ). SO2 is a known by-product of photochemistry7, which is generated when stellar ultra-violet radiation reacts with H2S and has been identified in WASP-107b at longer wavelengths with JWST/MIRI4. Finally, a narrow CH4 feature is identified at 3.32 μm (4.2σ). Although CH4 has been unobservable at both shorter and longer wavelengths2,4, the feature here is the Q-branch band head of CH4, which is strong enough to appear above the clouds and H2O (Fig. 2).

For H2O, we find volume-mixing ratio abundances of 10−1.85±0.22, which indicates super-solar abundances near 40× solar. The presence of CO2 is known to be a tracer of high metallicity15,16, and the retrieved abundances of 10−3.33±0.27 support this, although these values are lower than the chemical equilibrium expectations by about an order of magnitude as the molecule is further sensitive to atmospheric mixing. CH4 is also found to be heavily depleted, with retrieved abundances of 10−6.03±0.21. These CH4 abundances are consistent with previous Hubble Space Telescope (HST) and JWST upper limits2,4. At the millibar pressures probed in transmission spectra, at 40× solar the atmosphere in equilibrium would be about one part in a thousand CH4, whereas abundances slightly less than 1 part per million (ppm) are found. Measuring the amount of CH4 depletion is highly constraining for non-equilibrium models8, as it can be tied to the pressure and temperature of the hotter interior from which the species is mixed to higher levels. Although not every planet shows a CH4 depletion5, for those that do such as WASP-107b (ref. 4) and others3, only the upper limits on the depletion are available. By quantifying the CH4 depletion, the non-equilibrium chemistry can be tightly constrained and the depletion mechanisms can be explored. We find abundances of SO2 to be 10−5.06±0.13, which broadly match the JWST/MIRI results4 as well as those of WASP-39 b (refs. 7,17) indicating ongoing photochemistry in the atmosphere.

Given the retrieved molecular abundances, we ran a grid of forward atmospheric models to place constraints on the photochemistry, vertical mixing, metallicity and temperature structure of the planet. The model includes non-equilibrium chemistry from both vertical mixing and photochemistry self-consistently calculated18,19 (Methods). The chemistry was computed using T–P profiles in radiative–convective equilibrium, with a range of intrinsic temperatures, Tint. The intrinsic temperature is related to the emitted flux generated from the interior of the planet, which passes through the atmosphere, with Jupiter having a Tint near 100 K. With vertical mixing, most of the molecular species detected are expected to be uniformly mixed from their quench pressures at deeper and hotter conditions in chemical equilibrium to higher altitudes. We modelled the vertical mixing in a turbulent flow using the vertical eddy diffusion coefficient, Kzz. Our best-fit forward model to the retrieved abundances is shown in Fig. 3 (see also Extended Data Fig. 7). We find the combination of H2O, CO, CO2, CH4 and SO2 to be highly constraining. H2O and CO are insensitive to non-equilibrium chemistry in this temperature regime (Fig. 3) helping to uniquely measure the metallicity, whereas the combination of CO2 and CH4 are jointly sensitive to Kzz and Tint. SO2 also helps to provide an upper bound on Kzz (Extended Data Fig. 7). With these combined constraints, the forward model grid indicates the planet has a hot intrinsic temperature of Tint = 460 ± 40 K, a high metallicity of Z = 43 ± 8× solar and vigorous vertical mixing with a Kzz = 1011.6±0.1 cm2 s−1.

Fig. 3: Model interpretation.
figure 3

ATMO non-equilibrium chemistry model with vertical mixing and photochemistry. The best-fit non-equilibrium chemistry model (solid lines) and the abundance profiles in equilibrium (dot-dashed lines) are shown. The retrieved JWST abundances are shown (data points) with the model grid indicating that the planet has hot interior temperatures (Tint = 458 ± 38 K), a super-solar metallicity (Z/Z⊙ = 43 ± 8) with vigorous vertical mixing (Kzz = 1011.6±0.1 cm2 s−1).

Source Data

Full size image

Theoretical models have estimated the strength of vertical mixing and Kzz (refs. 20,21,22). However, the overall uncertainty on Kzz for giant planets remains large23, with values often ranging from 104 cm2  s−1 to 1012 cm2 s−1 for hot Jupiters8. The empirical constraints on Kzz for exoplanets currently rely on broadband Spitzer/IRAC photometry9, in which there is considerable planet-to-planet stochasticity and the molecular features are spectroscopically unresolved, giving rise to modelling degeneracies. The Kzz measurement in WASP-107b here is a substantial improvement over the order of magnitude estimates previously explored using either the upper limits on CH4 from JWST/MIRI4 or population-level results9. The measurement is one of the best evidence so far that non-equilibrium chemistry from vertical mixing is an important physical process in exoplanetary atmospheres, a mechanism long known to be important for Jupiter24. The Kzz value for WASP-107b at 1 bar is about 1,000–10,000 times higher than the typical estimates used7,25, which are largely based on three-dimensional general circulation models (GCM)20. This could indicate the planet has anomalously high vertical mixing, perhaps because of the high Tint resulting in a convective region that reaches lower than expected pressures. Alternatively, the GCMs may be inadequately capturing the eddy dissipation with the artificial viscosity parameters used26.

The high Tint is a direct constraint on the energy stored in the deep atmosphere, which can help us understand the mechanism responsible for the inflated radius of the planet. The anomalously large radii of irradiated warm Neptunes and hot Jupiters are one of the most intriguing and long-standing problems in our understanding of extrasolar giant planets with several proposed mechanisms, including tidal heating27, downward transport of kinetic energy28, enhanced opacities29, ongoing layered convection30, ohmic dissipation31 and, more recently, the advection of potential temperature32,33. The potential temperature mechanism has been recently demonstrated in the hot Jupiter WASP-76 b using first-principle radiative GCM simulations34. The hot interior temperature can explain the higher vertical mixing rates, as both hotter temperatures and lower surface gravities are expected to increase Kzz. These high vertical mixing rates with super-solar metallicities and hot intrinsic temperatures also explain the cloud properties observed, in which silicate cloud particles were found with JWST/MIRI4. For planets such as WASP-107b with equilibrium temperatures near 770 K, the condensation of silicate clouds should be occurring at high unobservable pressures if Tint was low (about 100 K) (ref. 8). However, with a high Tint of 460 K, the silicate cloud base is moved to about bar pressures (Extended Data Fig. 6) and the high Kzz values found here are more than sufficient to aloft the cloud particles to mbar pressures in which they are observed in transmission.

The mass, radius and Tint alone constrain the overall bulk metallicity of the planet to be \({Z}_{{\rm{p}}}={63.5}_{-8.5}^{+10.4} \% \). A precise atmospheric metallicity and Tint enable us to place better constraints on the overall metal content of the planet and estimate the core mass of the planet (Methods). So far, no gas giant exoplanet has had the presence of a core significantly detected, with upper limits reported for HAT-P-13 b (ref. 35) and WASP-107b (ref. 10). Any metals seen in the bulk but not in the atmosphere must be hidden in the interior in a core or composition layers. Assuming a uniform composition core, an isothermal 50:50% mixture of rock and water, gives us \({M}_{{\rm{c}}}={11.5}_{-3.6}^{+3.0}{M}_{\oplus }\)—which is one-third of the mass of the planet and the first statistically significant core detection for a giant exoplanet. This value is much higher than the previous estimates limiting the core to less than 4.6M⊕ (ref. 10), which had assumed the planet followed a standard cooling track with no tidal heating. These low core-mass values were in tension with standard core-accretion models11 that do not predict such a light core could accrete the massive H/He envelope observed. By contrast, our core-mass estimate matches the approximately 10M⊕ prediction needed to elicit runaway gas accretion within the lifetime of the disks. It is worth noting that the structure of the interior may not be exactly an envelope-on-core model, but may be more like Neptune and Uranus with a rocky core and an icy water mantle36 or the core may be diffuse and/or similar to that of Jupiter37,38. The core-mass result should be understood as the total excess non-H/He heavy elements in the interior, irrespective of how exactly it is structured. With the proof-of-concept demonstrated here, using CH4 depletion as a thermometer of the deep atmosphere for other gas giant exoplanets could help provide constraints at the population level about how planetary cores might differ with planet or host-star mass and metallicities.

WASP-107b is one of the lowest-density planets known, which has led to speculation that the planet has formed substantially different from the solar system, with dust-free accretion or in situ scenarios explored10. However, the important bulk properties found here line up well with the planets of the solar system (SS) indicating a similar core-accretion scenario. In particular, with a metallicity of 43 ± 8 × solar, WASP-107b is close to the mass-metallicity trend of the gas giant SS planets39. Moreover, with an estimated core mass of \({11.48}_{-3.6}^{+3.0}{M}_{\oplus }\), WASP-107b is also intermediate between Neptune and Saturn, falling reasonably along their lines of formation40,41. The main difference for WASP-107b seems to be its much hotter interior, resulting in an extremely puffy planet.

Methods

Data reduction

A transit of WASP-107b was observed on 23 June 2023 with JWST for 6.5 h using the NIRSpec instrument with the G395H grating. This setup gives a wavelength range from 2.7 μm to 5.18 μm with the corresponding resolving power ranging between 1,828 and 3,600. The instrument used the NRSRAPID readout pattern and SUB2048 subarray with 20 groups per integration, resulting in 1,230 integrations over the 6.5-h observation.

FIREFLy

We reduce the data using the Fast Infrared Exoplanet Fitting Lightcurve (FIREFLy)17,42 reduction suite, which starts with a customized reduction using the STScI pipeline and the uncalibrated images. Our custom reduction includes 1/f destriping at the group level before the ramp is fit. The jump-step and dark-current stages of the STScI pipeline were skipped. We then use the custom-run pipeline two-dimensional (2D) images after the gain scale step and perform customized cleaning of bad pixels, cosmic rays and hot pixels. Using cross-correlation, we measured the positional shift of the spectral trace across the detector and performed shift stabilization of the images with flux-conserving interpolation. This procedure has been found to reduce the amplitude of position-dependent trends17,42, as the underlying stellar spectra do not shift in wavelength across a pixel through the time series. We note that nearly identical results can be found if shift stabilization of the data is not performed and the extraction window is moved instead. Further position-dependent trends can be expected from intrapixel sensitivity variations13, although their magnitude on-sky has been very small as JWST generally has excellent pointing with the detector shifts found in this dataset to be of the order of about 1/1,000 of a pixel. We optimized the width of our flux-extraction aperture to minimize the standard deviation of the white-light curve photometry and extracted the spectrophotometric time series. The 2D extracted spectral time series can be seen in Extended Data Fig. 1. We extracted the transmission spectra with a range of different wavelength bin sizes, reporting an extraction with a resolution near R = 950.

We fit the transit light curves using a quadratic function to model stellar limb darkening given as

$$\frac{I(\mu )}{I(1)}=1-a(1-\mu )-b{(1-\mu )}^{2},$$
(1)

where I(1) is the intensity at the centre of the stellar disk, μ = cos(θ), where θ is the angle between the line of sight and the emergent intensity, and a and b are the limb-darkening coefficients. In practice, we remap a and b and fit instead for q1 and q2, which is shown to be less degenerate43. The resulting limb-darkening coefficients can be seen in Extended Data Fig. 3 with the best-fit orbital parameters given in Extended Data Table 1. Model limb-darkening coefficients calculated from three-dimensional models44 capture the wavelength dependence of limb darkening well, but we find an overall offset in the coefficients of −0.0378 and +0.175 for q1 and q2, respectively. We applied these offsets to the wavelength-dependent model limb-darkening coefficients and re-fit the transmission spectra with fixed limb darkening. This iterative procedure reduces the bin-to-bin scatter of the transmission spectra as the limb darkening is fixed, but allows for overall differences between the transit data and stellar model limb darkening. We bin the spectrophotometry into 576 wavelength bins, providing a spectral resolution of R ≈ 950 for the planet.

Stellar spot crossing

A stellar spot-crossing event is observed in our transit light curve data, appearing in both NRS1 and NRS2 shortly after mid-transit, which we decontaminated from the light curves and the transmission spectrum. The host star is known to be modestly active12, with optical photometric modulations of around 0.4% over a period of about 17 days. We modelled this spot empirically, first fitting the non-spotted region with a transit light curve and using the overall residual in the spotted region to define the photometric shape of the spot (Extended Data Fig. 2). A Gaussian filter was then applied to the spot shape and subsequently used to model the full transit light curves, fitting for a scaling factor. As the spot is low in amplitude and covers only a small portion of the data, our overall transmission spectrum is insensitive to the occulted spot. We found equivalent results between masking the region for all the transit light curve fits and fitting for the spot shape in the light curves.

TEATRO

An independent data reduction was performed using the Transiting Exoplanet Atmosphere Tool for Reduction of Observations (TEATRO) (https://github.com/ncrouzet/TEATRO) pipeline.

The data were processed using the jwst calibration software package v.1.10.2, CRDS v.11.16.19, CRDS context 1147 (refs. 45,46), starting from the ‘uncal’ files. For stage 1, a jump rejection threshold of 6 was used and the ‘jump.flag_4_neighbors’ parameter was turned off. For stage 2, the flat field and photometric calibration steps were skipped. Bad pixels and missing pixels were then corrected in each integration. The centre of the spectral trace in the spatial direction was computed by fitting a Gaussian function to each column and fitting a third-order polynomial to their maxima. For each integration, the background was computed on a column-by-column basis using pixels above and below the spectral trace (or only one of them at the short- and long-wavelength edges of NRS2). An aperture half-width of 2.73 and 3 pixels and background regions starting 6.5 and 6 pixels further away were used for nrs1 and nrs2, respectively. The background subtraction and spectral extraction were performed using the ‘extract_1d’ routine with the corresponding polynomial coefficients as parameters. The white-light flux was computed using the ‘white_light’ step routine. For the spectroscopic analysis, we binned the spectra in 86 and 136 wavelength bins of 0.01 μm width, covering the 2.86–3.72 μm and 3.82–5.18 μm ranges with nrs1 and nrs2, respectively. The light curves were normalized by the out-of-transit flux and a five-iteration 3.5 sigma-clipping was applied to remove the outliers.

Light-curve fits were performed using the exoplanet47,48 package with a quadratically limb-darkened transit light-curve model and a second-order polynomial to account for long-term trends. A fit to the white-light curve obtained from NRS1 only was first performed to derive system parameters (mid-transit time, impact parameter, planet-to-star radius ratio and stellar density). These parameters were then fixed for the spectroscopic light-curve fits, leaving only the radius ratio and polynomial trend as free parameters. Limb-darkening coefficients were obtained from ExoTETHyS49; they were free for the white-light curve fit and fixed for the spectroscopic light-curve fits. The median of the posterior distribution of the radius ratio was used to derive the transit depth in each wavelength bin. The uncertainties were computed by summing quadratically the standard deviations of the in- and out-of-transit parts of the light curve divided by the square root of their respective number of points (this gave a better estimate than the standard deviation of the posterior distributions).

Comparing the TEATRO reduction to FIREFLy, we find good consistency with the spectra agreeing at the point-to-point level at better than 1σ for both NRS1 and NRS2.

Resolution-linked bias

We have adjusted the retrieved abundance to account for the resolution-linked bias (RLB) effect50, in which planetary absorption signatures are diluted within stellar absorption lines. For G395H data, this, in particular, dilutes the transmission spectral signatures of CO, as the molecule is found in the atmosphere of both the star and planet. We followed the method in ref. 42, estimating the magnitude of the effect using high-resolution models. Compared with WASP-39A, WASP-107A is a cooler later-type star that enhances the RLB effect as the stellar CO is stronger. However, within the CO lines, the atmosphere of the planet is cloudier than WASP-39b, which reduces the effect. In total, we find that the RLB reduces the inferred CO abundance by 0.2 dex, which is about half of the uncertainty in the CO abundance.

Atmospheric models

ATMO setup

We use ATMO, a one-dimensional (1D)–2D radiative–convective equilibrium model for planetary atmospheres to generate models of WASP-107b. More comprehensive descriptions of the model can be found in refs. 18,32,51,52,53,54. The ATMO model is used here as both a forward physical model in radiative and chemical (dis-)equilibrium and a retrieval model in which the abundances and temperature–pressure (T–P) profile are let free to fit the data. By comparing the retrieved and forward model abundances within a single model suite, model-to-model systematics from such sources as differences in opacities can be avoided. ATMO has been validated against the publicly available MET office radiative transfer code SOCRATES55 and has been benchmarked against Exo-REM56 and petitCODE57 in the context of JWST observations58.

ATMO solves the radiative transfer equation with isotropic scattering in 1D plane-parallel geometry for the irradiation and thermal emission, finding a pressure–temperature (P–T) profile that satisfies hydrostatic equilibrium and conservation of energy and is also self-consistent with the atmospheric chemistry and opacities, given a set of elemental abundances. The total opacity of the gas mixture is computed using the correlated-k approximation using the random overlap method with resorting and rebinning51,59. The k-coefficients are calculated ‘on the fly’ for each atmospheric layer, spectral band and iteration such that the derived opacities are physically self-consistent with the T–P profile and chemical composition. ATMO can also be used as a full line-by-line code at high spectral resolution, although that is not used here as it is too computationally heavy for spectral retrieval purposes. The spectrally active species currently include H2–H2 and H2–He collision-induced absorption (CIA) opacities, as well as H2O, CO2, CO, CH4, NH3, Na, K, Li, Rb, Cs, TiO, VO, Fe, FeH, CrH, PH3, HCN, C2H2, H2S, SO2 and H– (see refs. 54,60 for full description). Multi-gas Rayleigh scattering contributions from the different species are also included.

The chemical abundances in equilibrium are determined by minimizing the Gibbs free energy following the method in ref. 61 and using the thermochemical data from ref. 62. Solar elemental abundances are set from refs. 63,64, with the chemistry fully flexible for any mix of input elemental abundances. This method also enables the depletion of gas phase species because of condensation as well as thermal ionization and dissociation. ATMO has several different chemical schemes that can be chosen and, by default, calculates the abundances for 175 neutral, 9 ionic and 93 condensate species. Condensation can be treated locally or with rainout. Rainout is treated by calculating the chemistry at the highest pressure initially, then following the T–P profile towards lower pressures. Condensed elements are removed locally, as well as for all lower pressures65, allowing for opacity changes that can alter the radiative–convective balance and P–T profile.

For non-equilibrium C–H–N–O chemistry, ATMO uses a chemical kinetics scheme from ref. 66 as described in ref. 18. As done in ref. 7 to model the S photochemistry of WASP-39b with ATMO, we used the thermochemical network from ref. 67 along with the photochemical scheme from ref. 68, including 71 photolysis reactions of H2S, S2, SO2, SO, SO2, CH3SH, SH, H2SO and COS.

ATMO includes a relatively simple treatment of clouds and hazes54 and does not consider the distribution of aerosol particles. An aerosol ‘haze’ scattering is implemented as enhanced Rayleigh-like scattering69, presented as

$$\sigma {(\lambda )}_{{\rm{haze}}}={\delta }_{{\rm{haze}}}{\sigma }_{{\rm{0}}}{(\lambda /{\lambda }_{0})}^{-{\alpha }_{{\rm{haze}}}},$$
(2)

where σ(λ) is the total scattering cross-section of the material, δhaze is an empirical enhancement factor, σ0 is the scattering cross-section of molecular hydrogen at 0.35 μm and αhaze is a factor determining the wavelength dependence with αhaze = 4 corresponding to Rayleigh scattering. Condensate ‘cloud’ absorption is assumed to have a grey wavelength dependence and is calculated as

$$\kappa {(\lambda )}_{{\rm{cloud}}}={\delta }_{{\rm{cloud}}}{\kappa }_{{{\rm{H}}}_{2}},$$
(3)

where κ(λ)cloud is the ‘cloud’ absorption opacity, δcloud is an empirical factor governing the strength of the grey scattering and \({\kappa }_{{{\rm{H}}}_{2}}\) is the scattering opacity due to H2 at 0.35 μm. σ(λ)haze and κ(λ)cloud are added to the total gaseous scattering and absorption, respectively.

For spectra retrievals, we coupled the forward ATMO model to a nested sampling statistical algorithm to marginalize the posterior distribution and measure the model evidence70,71,72. The retrieval aspects of ATMO were developed to fit transit and eclipse data, with results published for several transiting planets3,73,74,75,76,77. A main difference compared with the forward models is that the retrieval does not converge the T–P profile to radiative–convective equilibrium, but rather parameterizes the profile that is then fit against the planetary spectrum. With the retrieval model, the k-coefficients are still calculated ‘on the fly’ in the same way as described in the forward model (and equilibrium chemistry if specified) for every likelihood evaluation step to maximize accuracy and consistency. To parameterize the T–P profile, we use the three-channel (two optical, one infrared) analytic radiative equilibrium model as described in ref. 78, which is based on the derivations in ref. 79.

ATMO WASP-107b retrievals

For the NIRSpec WASP-107b data, we found the flexible parameterized P–T profile simply fit to a simple isotherm, so we instead fit for a single isothermal temperature for the atmosphere. Given the atmosphere of the planet is far out of equilibrium, we fit a constant VMR for each molecule. Apart from the cloud, we included the spectrally active species of H2–H2, H2–He (CIA) opacities, H2O (ref. 80), CO2 (ref. 81), CO (ref. 82), CH4 (ref. 83), NH3 (ref. 84), H2S (ref. 85) and SO2 (ref. 86).

We searched the G395H data for H2S and can marginally improve the model fits to the data, but the inclusion of the model is not supported with confidence by the Bayesian evidence. Considerable opacity from aerosols is also needed to fit the spectrum, which is in line with previous HST results in which H2O and He features are observed between 0.9 μm and 1.6 μm peaking above the clouds2,75. We fit the data with a grey cloud and a scattering haze but found only a grey cloud was needed to fit the data (Extended Data Table 3). Moreover, the choice of cloud parameterization did not have a strong influence on the retrieved abundances. Fitting the grey cloud either with a single uniform opacity or for a cloud-top pressure also gave similar results. A narrow feature of NH3 is potentially found in the data at 3.0 microns, but similar to H2S, its identification is not supported with high confidence by the Bayesian evidence. NH3 has been marginally detected at longer wavelengths with JWST/MIRI, indicating that the molecule could be confidently identified with a dedicated multi-wavelength study.

The retrieval constraints are given in Extended Data Table 2. CO is found to have abundances of 10−2.7±0.4, which is slightly low relative to 40× solar values derived from H2O. Although low, the uncertainty on the CO abundance is large as the feature is subtle in the data and H2O and clouds help mask the signature. A confounding factor regarding CO is an RLB effect17,50, where CO present in the stellar spectra dilutes the strong molecular line cores in the planetary transmission spectrum. This does not affect H2O, CH4 or CO2 as those molecules are not present in the stellar photosphere. In this work, we applied a correction for the CO abundance to account for the RLB effect.

Overall, we find a good fit to the data with a \({\chi }_{\nu }^{2}=1.1\) for 566 degrees of freedom and 10 free parameters. Our retrieval results are given in Extended Data Table 2 and Extended Data Fig. 4. To estimate the detection confidence of each molecule identified, we re-ran the retrieval leaving out the molecule in question and computed the detection significance using the Bayesian evidence between the models with and without the molecule. We report the results in which H2S and NH3 are included in the model, although the G395H data are not sufficient to detect either species confidently. We report the abundances and detection significances from the ATMO results, as they can be directly used to compare against the non-equilibrium models presented in the main text.

ATMO WASP-107b forward non-equilibrium chemistry models

It is computationally infeasible to currently run the ATMO retrieval with the full non-equilibrium photochemistry self-consistently computed at every model evaluation ‘on the fly’. Instead, we adopt a two-step grid-retrieval approach, first computing a grid of non-equilibrium chemistry values that are then fit against the retrieved VMR abundances. As we use the same ATMO model setup for both the forward and retrieval modelling, the abundances computed between them can be self-consistently compared.

Our use of a 1D atmospheric model to interpret the non-equilibrium chemistry has several assumptions. The first is that the atmosphere can be considered 1D and is independent of both latitude and longitude. Hot and ultra-hot Jupiters show large day–night temperature gradients87. Spitzer phase curves show these differences decrease at lower Teq (ref. 88). Thus, the latitude and longitude temperature differences for warm about 750 K planets such as WASP-107b are expected to be modest. Another important assumption is that the parameters of Tint and Kzz are constant and assumed to be uniform with depth. Our model contains a single convective region (Extended Data Fig. 6), although the possibility of a separate low-pressure convective zone has been studied89. Multiple convective regions would affect our interpretation of Tint and the link between atmospheric and interior constraints. Furthermore, our modelling also assumes a single Kzz value that is constant with altitude. Theoretical studies have predicted that Kzz increases with altitude20. Higher Kzz values at higher altitudes would not greatly affect our results, as the quench pressures for the molecular features probed here are similar (Extended Data Fig. 6). However, if multiple convective zones are present, the use of a single Kzz value may have to be re-visited.

To estimate the parameter space needed for the grid, we first used a chemical equilibrium model90 and assumed a single quench pressure to estimate the posterior distribution of temperature, pressure, metallicity and C/O ratio parameters that are consistent with the retrieved abundances (Extended Data Fig. 5). This equilibrium model indicated metallicities near 50× solar, C/O ratios near 0.5 and temperatures near 1,000 K at pressures just below a bar are needed to fit the observed molecular abundances. These temperatures are reached at 1 bar only if the planet has an intrinsic temperature greater than about 300 K with super-solar metallicities (Extended Data Fig. 6). Intrinsic temperatures lower than 300 K are ruled out as CH4 becomes too abundant in the deep atmosphere, such that no amount of mixing could sufficiently deplete the molecule. We also estimate a limit on Kzz to be less than 1013.5 cm2 s−1 based on requiring velocities to be less than the local sound speed. The presence of clouds near millibar pressures also places a lower limit on Kzz of about 109 cm2 s−1, as significant vertical mixing is needed to keep the aerosol particles aloft20.

We input P–T profiles in radiative–convective equilibrium, using intrinsic temperatures (Tint) of 300–500 K in steps of 50 K along with a Tint of 425 K. With these T–P profiles, we computed the non-equilibrium chemistry as described in ref. 18, using metallicities (Z) between 20×  solar and 50× solar in steps of 5×, and 10 log[Kzz (cm2 s−1)] values ranging from 11 to 13. A K6V star was used for our input to the photochemical model91. Each non-equilibrium model was integrated to 1 × 108 s, which we found was sufficient such that the abundances did not evolve further. We assumed a solar C/O ratio, which is also consistent with WASP-107A from ref. 92 which found values of C/O = 0.5 ± 0.1. We find the C/O ratio does not differ substantially from the solar or host-star value (Extended Data Fig. 5), and we subsequently ran low (0.1 and 0.2) and high (0.7) non-equilibrium chemistry cases that did not improve the fits.

For this grid of models, we then calculated the χ2 for each grid model, using the retrieved abundances for H2O, CO, CO2, CH4 and SO2 in Extended Data Table 2. We did not use the H2S or NH3 abundances, as they are not fully supported by the data and the errors in any case are large. The best-fit model has a good χ2 of 4.1 when fitting the five abundances with three grid parameters (Fig. 3 and Extended Data Fig. 7). The grid-posterior can be seen in Extended Data Fig. 7, in which each model was given a weight proportional to its probability, P, calculated as \(P={{\rm{e}}}^{-{\chi }^{2}/2}\). From the marginalized grid-posterior, we fit a Gaussian to the distributions and report the fitted Z/Z⊙, Tint and Kzz values in Extended Data Table 2.

We find the molecular constraints highly constraining for both Tint and Kzz. Lower vertical mixing rates, generally, require higher intrinsic temperatures to deplete CH4 to similar values. However, if Tint becomes too high, it over-depletes CH4 relative to CO2. Moreover, SO2 helps provide an upper bound on Kzz as SO2 is sensitive to both the photochemistry, needed to produce the species at low pressures, and vertical mixing. Higher values of Kzz mix the species from deeper pressures at which the abundances are lower, whereas lower Kzz values allow the photochemistry to build up the species near 0.1 mbar pressures.

NEMESIS free retrievals

NEMESIS is a free-chemistry radiative transfer and retrieval code, initially developed for solar system applications93. It couples a correlated-k59 radiative transfer scheme with the PyMultiNest70,71,72,94 Nested Sampling algorithm95. It has been extensively applied to spectra of hot Jupiter exoplanets96,97,98.

The retrievals for WASP-107b incorporate abundances of the gases H2O (line data99), CO2 (ref. 100), CO (ref. 101), CH4 (ref. 102), SO2 (ref. 86), NH3 (ref. 103), HCN (ref. 104) and H2S (ref. 105). For each of the gas abundances, we specify a log-uniform prior spanning 10−12–10−0.5. The rest of the atmosphere is composed of H2 and He with a ratio of 0.8547:0.1453. CIA for H2 and He is taken from refs. 106,107.

Other retrieved properties are an isothermal temperature with a uniform prior and a range of 300–1,000 K, and a grey cloud-top pressure with a log-uniform prior and a range of 10−8–100 bars. We also retrieve the planetary radius at a reference pressure of 100 bars. The retrieved abundances are given in Extended Data Table 2. Additional retrievals were run that had a more flexible T–P profile along with retrievals with patchy clouds and enstatite clouds, although they were not favoured statistically over the simpler model.

Interior structure models

The constraints we found for the intrinsic temperature of WASP-107b have important implications for the interior structure of the planet. To investigate these, we use the interior structure models of ref. 39, revised to use the H/He EOS108 and parameterize the amount of metal in the core versus mixed into the envelope. As in ref. 39, we assume the metal is a 50–50 mixture of rock and ice.

To match our models to WASP-107b, we use a Bayesian retrieval approach similar to the one used in ref. 109. However, our measurement of Tint determines the thermal state of the plane rather than thermal evolution. The parameters of the statistical model were, therefore, the mass of the planet M, the metallicity of the envelope Ze (a mass fraction), the core mass Mc and the specific entropy of the envelope s. From these, our forward model calculates the radius and Tint of the planet. The model was constrained by four normally distributed observations: the mass, the radius, the atmospheric metallicity and Tint. These are constraining enough that we can choose very uninformative priors to keep the parameters inside model bounds: a uniform mass prior from 0.01 to 30 MJ, a uniform Zenv from 0 to 1, a uniform core mass from 0 to MJ (thus conditional on M but still proper) and a uniform entropy prior from 5 to 11 kb per baryon. We sampled from the posterior using the Metropolis–Hastings Markov chain Monte Carlo algorithm and verified convergence with the auto-correlation plots and the Gelman–Rubin statistic110.

Our interior model retrievals (Extended Data Fig. 8) find a planet that contains similar amounts of H/He compared with heavier elements, with Zp = 0.608 ± 0.046. This seems reasonable for a planet of this mass: it accreted a substantial amount of material from the disk, but never reached the runaway gas accretion stage. However, this is in contrast to the previous estimates of the composition of the planet, which saw it as having substantially larger quantities of H/He: that is, more than 85% H/He from ref. 10. This is because our measured Tint is much hotter than a standard thermal evolution model would suggest; we would expect the planet to have cooled below our 2σ lower bound of 300 K within tens of Myr of formation. Instead, it is clear that the planet is being inflated by an extra heat source despite being too cool to experience hot Jupiter-type inflation111. Instead, a mechanism such as tidal heating112,113,114, previously suspected to be operating on WASP-107b (refs. 10,115), is warming the interior. These models propose that the eccentricity of WASP-107b is excited by WASP-107 c, but then damped out through tidal interaction with the star; the resulting energy is deposited in the interior of the planet. Our large measured Tint is evidence that this process could be occurring, with a precise eccentricity further capable of constraining the mechanism.

Another interesting result of our interior structure models is that it puts constraints on the mass of the core of the planet. The mass, radius and Tint already constrain the overall metallicity of the planet (more metal means a smaller planet) to \({Z}_{{\rm{p}}}=0.63{5}_{-0.085}^{+0.104}\), but we also have a measurement of the atmospheric metallicity of the planet. This enables us to make a modestly better constraint on Z (see previous paragraph); more importantly, metal seen in the bulk but not in the atmosphere must be hidden in the interior in a core or composition layers (see ref. 109 for further discussion). Assuming a uniform composition core gives us \({M}_{{\rm{c}}}={11.5}_{-3.58}^{+3.01}{M}_{\oplus }\). This range excludes zero, meaning that we have detected the presence of a core. The exact structure of the interior may not be exactly an envelope-on-core model, but may be more like Neptune and Uranus with a rocky core, water envelope, than a H/He envelope36. Alternatively, the core may be diffuse and/or layered as the core of Jupiter is thought to be37,38. Nevertheless, our result that a core exists is not affected by these possibilities: we cannot propose a layered core without a core in the first place.