# A high internal heat flux and large core in a warm Neptune exoplanet

Luis Welbanks, Taylor J. Bell, Thomas G. Beatty, Michael R. Line, Kazumasa Ohno, Jonathan J. Fortney, Everett Schlawin, Thomas P. Greene, Emily Rauscher, Peter Mcgill, et al.

► To cite this version:

Luis Welbanks, Taylor J. Bell, Thomas G. Beatty, Michael R. Line, Kazumasa Ohno, et al. A high internal heat flux and large core in a warm Neptune exoplanet. Nature, 2024, 630, pp.836-840. 10.1038/s41586-024-07514-w. insu-04843530

HAL Id: insu-04843530

https://insu.hal.science/insu-04843530v1

Submitted on 4 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A high internal heat flux and large core in a warm Neptune exoplanet

https://doi.org/10.1038/s41586-024-07514-w

Received: 15 January 2024

Accepted: 2 May 2024

Published online: 20 May 2024

<figure><img src="image_2.png" /></figure>

Check for updates

Luis Welbanks¹²,³, Taylor J. Bell²,³, Thomas G. Beatty⁴, Michael R. Line<fnref n="1" />, Kazumasa Ohno⁵,⁶, Jonathan J. Fortney⁵, Everett Schlawin⁷, Thomas P. Greene³, Emily Rauscher⁸, Peter McGill⁹, Matthew Murphy⁷, Vivien Parmentier¹⁰, Yao Tang⁵, Isaac Edelman², Sagnick Mukherjee⁵, Lindsey S. Wiser<fnref n="1" />, Pierre-Olivier Lagage¹¹, Achrène Dyrek¹¹ & Kenneth E. Arnold⁴

Interactions between exoplanetary atmospheres and internal properties have long been proposed to be drivers of the inflation mechanisms of gaseous planets and apparent atmospheric chemical disequilibrium conditions<fnref n="1" />. However, transmission spectra of exoplanets have been limited in their ability to observationally confirm these theories owing to the limited wavelength coverage of the Hubble Space Telescope (HST) and inferences of single molecules, mostly H₂O (ref. 2). In this work, we present the panchromatic transmission spectrum of the approximately 750 K, low-density, Neptune-sized exoplanet WASP-107b using a combination of HST Wide Field Camera 3 (WFC3) and JWST Near-Infrared Camera (NIRCam) and Mid-Infrared Instrument (MIRI). From this spectrum, we detect spectroscopic features resulting from H₂O (21σ), CH₄ (5σ), CO (7σ), CO₂ (29σ), SO₂ (9σ) and NH₃ (6σ). The presence of these molecules enables constraints on the atmospheric metal enrichment (M/H is 10⁻¹⁸× solar³), vertical mixing strength (log₁₀ Kzz = 8.4–9.0 cm² s⁻<fnref n="1" />) and internal temperature (≥345 K). The high internal temperature is suggestive of tidally driven inflation⁴ acting on a Neptune-like internal structure, which can naturally explain the large radius and low density of the planet. These findings suggest that eccentricity-driven tidal heating is a critical process governing atmospheric chemistry and interior-structure inferences for most of the cool (<1,000 K) super-Earth-to-Saturn-mass exoplanet population.

The mass of WASP-107b is similar to Neptune (1.78 M<sub>⊕</sub>), but its extreme low density is suggestive of a hydrogen/helium (H/He) envelope-to-core-mass ratio (>85% (ref. 5)) more like Jupiter/Saturn (around 90%) than like Neptune/Uranus (5–15% (ref. 6)). This high an envelope mass fraction for such a low-mass planet presents challenges to the standard core-accretion model of planet formation⁷—it is unclear how such a low-mass planetary core (inferred to be <5 M<sub>⊕</sub> (ref. 5)) could accrete such a massive gaseous envelope, but then stop short of fully growing into a 'Jupiter'⁸. An alternative hypothesis⁴ is that tidal heating could inflate the planetary envelope, reducing the need for such an envelope-to-core-mass ratio.

As well as mass and radius constraints, measurements of the atmospheric abundances of molecules containing carbon (C), oxygen (O), nitrogen (N) and sulfur (S) can be used to jointly constrain both scenarios. The abundances of molecules in the atmosphere will be set by the intrinsic elemental inventory and the chemical processes primarily driven to disequilibrium induced by transport and photochemistry⁸,⁹. The former constrains the partitioning of material between the envelope and core, whereas the latter is sensitive to the amount of tidal heating, altering the deep-atmosphere temperature and, consequently,

the molecular abundances at their quenched values<fnref n="1" />. One of the key challenges in interpreting exoplanet atmosphere compositions is disentangling the intrinsic elemental inventory from different chemical processes, given the presence (or lack thereof) of molecular spectral features.

Initial reconnaissance spectroscopy with the HST/WFC3 (0.8–1.6 µm) presented²,¹⁰ a muted water vapour absorption feature (1.4 µm) and a surprising lack of methane (CH₄) absorption (1.6 µm)—expected to be in abundance under solar elemental abundances and chemical equilibrium at temperatures below approximately 1,000 K (refs. 9, 11). The presence of water but lack of methane could either indicate a low carbon-to-oxygen ratio (C/O) envelope or be because of the quenching of methane at deeper, hotter layers<fnref n="1" />,². The constraints on the water abundance were not precise enough to determine the bulk envelope metal enrichment.

The broad wavelength coverage (0.6–28 µm) offered by the JWST provides access to several molecular bands of the major C, O, N and S species, critical to enabling the precise atmospheric abundance constraints necessary for breaking the aforementioned degeneracies. As part of the MANATEE NIRCam + MIRI Guaranteed Time Observations

<footnote n="1">¹School of Earth and Space Exploration, Arizona State University, Tempe, AZ, USA. ²Bay Area Environmental Research Institute, NASA's Ames Research Center, Moffett Field, CA, USA. ³Space Science and Astrobiology Division, NASA's Ames Research Center, Moffett Field, CA, USA. ⁴Department of Astronomy, University of Wisconsin–Madison, Madison, WI, USA. ⁵Department of Astronomy and Astrophysics, University of California, Santa Cruz, Santa Cruz, CA, USA. ⁶Division of Science, National Astronomical Observatory of Japan (NAOJ), Tokyo, Japan. ⁷Steward Observatory, University of Arizona, Tucson, AZ, USA. ⁸Department of Astronomy, University of Michigan, Ann Arbor, MI, USA. ⁹Space Science Institute, Lawrence Livermore National Laboratory, Livermore, CA, USA. ¹⁰Laboratoire Lagrange, Observatoire de la Côte d'Azur, Université Côte d'Azur, Nice, France. ¹¹Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, Gif-sur-Yvette, France. ¹²e-mail: luis.welbanks@asu.edu</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1 | Spectroscopic and broadband NIRCam light curves of the transit of WASP-107b.** The raw spectroscopic F322W2 and F444W transit light curves are shown in panels a and b after spectral binning (0.015-µm bins) but without any temporal binning. Masked values (for example, cosmic rays) have been coloured black. Even from these raw data, it is possible to visually identify the increased transit depth ≲3 µm primarily caused by H₂O, from 3.9-4.1 µm caused by SO₂.</figcaption></figure>

and from 4.3-4.5 µm caused by CO₂. The F322W2 broadband (2.45-3.95 µm) and F444W broadband (3.89-4.97 µm) transit light curves are shown in panels c and d with grey points without error bars. Black points with 1σ error bars show temporally binned data with a cadence of 5 min. Note that the error bars are typically smaller than the point size. BJD<sub>TDB</sub> is the date in the Barycentric Julian Date in the Barycentric Dynamical Time system.

(GTO) programme (JWST-GTO-1185 (ref. 12)), we collected two new transit observations of WASP-107b using JWST/NIRCam's F322W2 (2.4-4.0 µm) and F444W (3.9-5.0 µm) filters and grism¹³ on 14 January 2023 and 4 July 2023, respectively. We analysed the JWST/NIRCam observations with three separate data-analysis pipelines (Eureka!¹⁴, Pegasus (Beatty et al., in preparation) and tshirt¹⁵), which all agree well within error (Fig. 1 and Extended Data Fig. 2); we ultimately chose to adopt the analysis of Eureka! for our fiducial NIRCam spectrum when performing atmospheric modelling. Also, we incorporate recently published JWST/MIRI Low-Resolution Spectrometer

(LRS) observations¹⁶, as well as the previously published HST/WFC3 G102 and G141 observations²,¹⁰. To ensure a uniform set of orbital and limb-darkening parameters between the many different instruments, we re-analysed the HST/WFC3 observations using Pegasus and the JWST/MIRI LRS observations using Eureka! We find that these re-analysed WFC3 and MIRI/LRS spectra are consistent with the literature spectra except for a constant offset caused by the improved orbital solution. Overall, these data give us a panchromatic dataset with observations spanning 0.8-12.2 µm, with continuous coverage between 2.45 and 12.2 µm.

<figure><img src="image_4.png" /><figcaption>**Fig. 2 | Independent reductions of the transmission spectrum of WASP-107b.**
*a*, The NIRCam F322W2 and F444W transmission spectra as reduced by the Eureka!, Pegasus and tshirt pipelines are shown in different colours at a constant wavelength binning of $\Delta\lambda = 0.015$ µm and with 1$\sigma$ error bars. All three analyses show clear agreement on the overall shape of the spectrum, with features from H₂O (approximately 2.5-3.5 µm), CH₄ (approximately 3.2-3.8 µm), SO₂ (approximately 3.9-4.1 µm) and CO₂ (approximately 4.2-4.6 µm) all clearly visible by eye in each of the reductions. **b**, The absorption cross-sections of the six detected species, many of which are visually identifiable in panel **a**.</figcaption></figure>
<figure><img src="image_5.png" /></figure>

**Fig. 3** | Interpretation of the transmission spectrum of WASP-107b. The observed transmission spectrum of the planet with $1\sigma$ error bars is compared with the best-fit 1D-RCPE model ($\chi^2/N_{data} = 1.4$) shown at a spectral resolution of $R = 300$. The coloured shaded regions show the contributions from individual gases to the best-fit model. The dashed grey line shows the clear-atmosphere component of the model, that is, the gas contributions without the presence of

clouds or hazes. Panel a shows the broad-wavelength spectrum (0.8–12 µm) of the planet. Panel b shows the fit residuals with the corresponding $2\sigma$ data error envelope. Panel c shows the NIRCam observations, whereas panel d shows HST/WFC3 spectra and panel e shows JWST/MIRI observations. The HST/WFC3 data are binned at a constant $\Delta\lambda = 0.025$ µm, the NIRCam data are binned at $\Delta\lambda = 0.015$ µm and the MIRI data are binned at $\Delta\lambda = 0.15025$ µm.

Among the chemical species expected in exoplanet atmospheres⁹¹⁷ (see Methods), the final NIRCam spectrum shown in Fig. 2 shows prominent absorption features owing to H₂O (2.5–3.2 µm, detected at 21σ), CO₂ (2.66–2.86 µm and 4.2–4.5 µm, detected at 29σ), CO (4.5–4.9 µm, detected at 7σ) and SO₂ (3.94–4.1 µm, detected at 9σ), with weaker features owing to CH₄ (3.2–3.5 µm, detected at 5σ) and NH₃ (2.9–3.1 µm, detected at 6σ), with reported detections from the one-dimensional radiative-convective-photochemical equilibrium (1D-RCPE) models described below. These features are complemented by two further H₂O bands in HST/WFC3 (around 1.13 µm and 1.4 µm) and another in MIRI (between 5 and 7 µm), along with another strong SO₂ feature (7.1–7.7 µm) and a weak NH₃ feature (10.3–11 µm) in MIRI, shown in Fig. 3. The panchromatic spectrum exhibits a slight downward slope from blue to red, indicative of aerosol scattering¹⁸⁻²⁰, and a strong concavity across MIRI, previously attributed to silicate cloud particulate resonance features¹⁶. The presence of both CO₂ and SO₂ features is indicative of envelope metal enrichment above solar and photochemistry²¹, whereas the relatively weak CH₄ feature is suggestive of CH₄ depletion²,¹⁰,¹⁶.

To rigorously infer the atmospheric composition and internal temperature from the observed transmission spectrum, we use Bayesian inference with the transmission spectra generated from a suite of

1D-RCPE models¹¹ (see Methods). This self-consistent method correctly captures the degeneracies between the intrinsic atmospheric composition and chemical processes. The final result from this process are samples from a posterior probability distribution that represents the constraints on $T_{itr}$, $T_{int}$, [M/H], C/O and cloud properties, and an estimate of the Bayesian model evidence, which we use to draw conclusions about the atmosphere of WASP-107b. Figure 3 shows that this model setup adequately explains the observed panchromatic spectrum, capturing nearly all of the salient features. From this process, we find a metallicity of 10–18× solar³ (at 68% confidence, median of 12.3× solar) and a sub-solar carbon-to-oxygen ratio (C/O = 0.33±0.05). Furthermore, the spectrum also requires a high internal temperature ($T_{int} > 345$ K at 99.7% confidence, whereas a value of <100 K would be expected⁴ given the low mass of the planet and the approximately 3-Gyr age of the star³) and an atmosphere with a strong enough vertical mixing (that is, eddy diffusion, $\log K_{zz} = 8.6^{+0.4}_{-0.2}$) to quench methane along the deeper (0.25–0.65 bar) and hotter (around 1,100–1,300 K) parts of the atmosphere, confirming previous suggestions²,¹⁶.

We further validate the interpretation from the self-consistent 1D-RCPE models by performing more Bayesian inferences using parametric atmospheric models that independently fit for the chemical
<figure><img src="image_6.png" /><figcaption>**Fig. 4** | Inferred molecular volume mixing ratios from the transmission spectrum of WASP-107b. The retrieved volume mixing ratios (that is, abundances) for each detected gas (panels **a–f**) are shown from the free retrieval (grey posterior distributions, with median and 68% confidence error bars) and from the 1D-RCPE grid inference (median and 68% confidence profiles shown as coloured pressure-dependent abundance profiles). The solid black lines show expectations for an atmosphere in chemical equilibrium under the inferred [M/H] and C/O. The dotted black lines show expectations for chemical disequilibrium under a cool internal temperature of 200 K. The self-consistent abundance profiles from the 1D-RCPE retrieval for a high internal temperature are generally consistent with the retrieved abundances from the free retrieval within 1σ. The shaded red region spans pressures from 1 mbar to 10⁻⁵ bar and corresponds to the pressures examined by our observations. The colours of the abundance profiles match the contributions in Fig. 3 and the cross-sections in panel **b** of Fig. 2.</figcaption></figure>

abundances and vertical pressure–temperature structure of the planetary atmosphere without any assumptions of radiative–convective thermo-chemical equilibrium (called ‘free retrievals’). Using two independent inference frameworks (Aurora and CHIMERA; see Methods, detection significances from CHIMERA), we confirm the detections of CO₂ (27σ), H₂O (18σ), SO₂ (8σ), CO (5σ), NH₃ (5σ) and CH₄ (8σ) at strong confidence, in agreement with the 1D-RCPE models. We find that the detection of NH₃ is mostly driven by the NIRCam observations and do not depend only on the red edge of the MIRI observations (see Methods). The derived chemical abundances, assumed to be constant with height by these models, are also consistent with the inferred abundance profiles from the 1D-RCPE models (see Fig. 4) and are also suggestive of a metal-enriched envelope with depleted CH₄ abundances.

The relatively high internal temperature we infer for WASP-107b is probably caused by tidal heating in the planetary interior. Recent radial-velocity observations measured the planet to have a mildly non-zero orbital eccentricity of *e* = 0.06 ± 0.04 (ref. 5). Tidal-heating relations²² predict that, at this nominal *e* = 0.06, the internal temperature of WASP-107b will be *T*<sub>int</sub> ≈ 350 K for a Neptune-like tidal quality factor of *Q* = 10⁴. The expected internal temperature from tidal heating decreases to *T*<sub>int</sub> ≈ 200 K to *T*<sub>int</sub> ≈ 110 K for Jupiter-like tidal quality factors of *Q* = 10⁵ to *Q* = 10⁶. Further heating of the deep atmosphere may be achieved by vertical and horizontal mixing²³⁻²⁶. If we assume that tidal effects are solely responsible for heating the atmosphere of WASP-107b, then our retrieved 3σ lower limit of *T*<sub>int</sub> > 345 K implies a corresponding 3σ upper limit on the tidal quality factor of the planet of *Q* < 10³.⁸. If WASP-107b has a tidal quality factor much lower than Neptune (10³.⁹ ≲ *Q*<sub>N</sub> ≲ 10⁴.⁵)²⁷, this would be consistent with the inference of a large core for the planet, as typical core material for ice giants is expected to have *Q* ≈ 10², whereas typical gas envelopes have *Q* ≈ 10⁵ (ref. 28).

The large value of *T*<sub>int</sub> implies a much hotter, lower-density H/He envelope that has been previously appreciated. Using state-of-the-art structure models²⁹, along with *T*<sub>int</sub> = 350 K (the retrieved 3σ lower limit), we find that the radius of the planet can be explained by a model that has >22 *M*<sub>⊕</sub> of rock/iron in its interior, modelled here as a distinct core. This ratio of solids to H/He is similar to that of the bulk composition of Uranus and Neptune in our Solar System³⁰. Our revised view of the planet is different from previous work⁵,³¹, as without tidal heating, the

low bulk density of the planet can only be explained by a structure that is mostly H/He with a very small core (<5 *M*<sub>⊕</sub> (ref. 5)), perhaps at odds with core-accretion theory.

The inferred atmospheric metallicity of WASP-107b is lower than expectations from the Solar System metal-enrichment trend³² as determined by the CH₄ abundances³³, which predicts an enhancement of about 32× solar (see also ref. 34). Our results confirm that different elements in a planetary atmosphere can be differently enhanced, as previously suggested by HST observations³². Furthermore, our results demonstrate that individual molecules (for example, H₂O) may not be good bulk metallicity tracers. Instead, JWST broadband spectra, as analysed here, present a key opportunity to derive estimates of the bulk atmospheric metal enrichment informed by the measured abundances of several gases and by self-consistently considering the impact of disequilibrium processes arising from interactions between the interior and the atmosphere of the planet.

The understanding of vertical mixing in the atmospheres of giant planets and brown dwarfs remains a notable challenge in atmospheric physics¹,³⁵. The inferred vertical mixing strength (*K*<sub>zz</sub>) can be high (*K*<sub>zz</sub> ≈ 10<sup>8</sup>–10<sup>11</sup> cm<sup>2</sup> s<sup>-1</sup>, depending on *T*<sub>int</sub>) in convective zones in which mixing is driven by convective overturn and mixing-length theory can be used, or much more uncertain in radiative zones in which the driving mechanism remains unclear³⁶,³⁷. Recent studies of mixing in brown-dwarf atmospheres, at temperatures similar to WASP-107b (refs. 35,38), infer low *K*<sub>zz</sub> values of 10<sup>2</sup>–10<sup>5</sup> cm<sup>2</sup> s<sup>-1</sup> owing to sluggish mixing in deep atmospheric radiative zones in which chemical abundances are quenched. Studies of this same phenomenon in Jupiter have long suggested a *K*<sub>zz</sub> value of approximately 10<sup>8</sup>–10<sup>9</sup> cm<sup>2</sup> s<sup>-1</sup>, although recent work³⁹ points to a lower value (*K*<sub>zz</sub> ≈ 10<sup>6</sup> cm<sup>2</sup> s<sup>-1</sup>) as the best explanation for the CO and H₂O abundances of Jupiter, again perhaps because of a deep radiative zone. Overall, the high *K*<sub>zz</sub> values implied for WASP-107b are notable in comparison with other objects. This suggests the quenching of chemical abundances in an atmospheric convective zone, with a high internal flux and shallow radiative–convective boundary. This can only be achieved if substantial extra internal energy, such as from tidal heating, keeps the interior and deep atmosphere much warmer than a standard cooling model would suggest.

The detection of the major C, O, N and S reservoirs in the spectrum of WASP-107b resulting from the broad wavelength coverage
demonstrates the unparalleled capabilities of the JWST for the detailed atmospheric characterization of exoplanet atmospheres. The simultaneous constraints on several chemical species enables unique solutions for the bulk atmospheric metal enrichment and strength of vertical mixing and provides strong evidence for a high internal temperature most probably arising from tidal heating owing to a small, yet notable, orbital eccentricity. The combination of constraints on atmospheric metallicity and the internal heat flux together provide new constraints on the relative gas-to-core mass fraction, naturally explaining the low planetary density. These initial constraints on the elemental ratios of C, O, S and N have the potential to inform the planetary accretion history⁴⁰.

The detection and constraint on the abundance of CH₄ in this warm Saturn adds to the growing number of inferences of this sought-after carbon-bearing species in transiting exoplanets¹¹,⁴¹. Furthermore, the low abundance of methane relative to thermochemical expectations confirms its sensitivity to mixing processes and the internal temperature/heat flux¹, providing a critical atmospheric diagnostic of interior processes and structure. We anticipate that the most common type of planet across the galaxy will have elevated internal temperatures arising from tidal heating and thus should have relative depletions of methane. Of the 262 known cool (<1,000 K) 'Neptune-like'⁴² planets (with well-determined mass and radii), approximately two-thirds have a reported non-zero eccentricity⁴³. As with WASP-107b, deciphering the intrinsic nature of these worlds will be critically dependent on the constraints on the total elemental inventory, the strength of vertical mixing and the internal temperature. Future studies from the MANATEE programme will explore the transmission and emission spectrum of warm exoplanets to further disentangle the effects of disequilibrium chemistry owing to interior-atmosphere interactions, photochemistry and directly link this population to the Solar System gas giants and their theorized formation pathways.

## Online content

Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-024-07514-w.

1. Fortney, J. J. et al. Beyond equilibrium temperature: how the atmosphere/interior connection affects the onset of methane, ammonia, and clouds in warm transiting giant planets. *Astron. J.* **160**, 288 (2020).

2. Kreidberg, L., Line, M. R., Thorngren, D., Morley, C. V. & Stevenson, K. B. Water, high-altitude condensates, and possible methane depletion in the atmosphere of the warm super-Neptune WASP-107b. *Astrophys. J. Lett.* **858**, L6 (2018).

3. Lodders, K., Palme, H. & Gail, H. P. in Landolt-Börnstein - Group VI Astronomy and Astrophysics, Vol. 4B (ed. Trümper, J. E.) 712 (Springer, 2009).

4. Millholland, S., Petigura, E. & Batygin, K. Tidal inflation reconciles low-density sub-Saturns with core accretion. *Astrophys. J.* **897**, 7 (2020).

5. Piaulet, C. et al. WASP-107b's density is even lower: a case study for the physics of planetary gas envelope accretion and orbital migration. *Astron. J.* **161**, 70 (2021).

6. Guillot, T. et al. in Protostars and Planets VII (eds Nutsuka, S., Aikawa, Y., Muto, T., Tomida, K. & Tamura, M.) 947 (Astronomical Society of the Pacific, 2023).

7. Pollack, J. B. et al. Formation of the giant planets by concurrent accretion of solids and gas. *Icarus* **124**, 62–85 (1996).

8. Madhusudhan, N., Knutson, H., Fortney, J. J. & Barman, T. in Protostars and Planets VI (eds Beuther, H., Klessen, R. S., Dullemond, C. P. & Henning, T.) 739 (Univ. Arizona Press, 2014).

9. Moses, J. I. et al. Compositional diversity in the atmospheres of hot Neptunes, with application to GJ 436b. *Astrophys. J.* **777**, 34 (2013).

10. Spake, J. J. et al. Helium in the eroding atmosphere of an exoplanet. *Nature* **557**, 68–70 (2018).

11. Bell, T. J. et al. Methane throughout the atmosphere of the warm exoplanet WASP-80b. *Nature* **623**, 709–712 (2023).

12. Schlawin, E., Greene, T. P., Line, M., Fortney, J. J. & Rieke, M. Clear and cloudy exoplanet forecasts for JWST: maps, retrieved composition, and constraints on formation with MIRI and NIRCam. *Astron. J.* **156**, 40 (2018).

13. Greene, T. P. et al. $\lambda = 2.4$ to 5 µm spectroscopy with the James Webb Space Telescope NIRCam instrument. *J. Astron. Telesc. Instrum. Syst.* **3**, 035001 (2017).

14. Bell, T. et al. Eureka!: an end-to-end pipeline for JWST time-series observations. *J. Open Source Softw.* **7**, 4503 (2022).

15. Schlawin, E. & Glidic, K. tshirt. GitHub https://github.com/eas342/tshirt (2022).

16. Dyrek, A. et al. $SO_2$, silicate clouds, but no CH₄ detected in a warm Neptune. *Nature* **625**, 51–54 (2024).

17. Madhusudhan, N. Exoplanetary atmospheres: key insights, challenges, and prospects. *Annu. Rev. Astron. Astrophys.* **57**, 617–663 (2019).

18. Lecavelier Des Etangs, A., Pont, F., Vidal-Madjar, A. & Sing, D. Rayleigh scattering in the transit spectrum of HD 189733b. *Astron. Astrophys.* **481**, L83–L86 (2008).

19. Sing, D. K. et al. A continuum from clear to cloudy hot-Jupiter exoplanets without primordial water depletion. *Nature* **529**, 59–62 (2016).

20. Ohno, K. & Kawashima, Y. Super-Rayleigh slopes in transmission spectra of exoplanets generated by photochemical haze. *Astrophys. J. Lett.* **895**, L47 (2020).

21. Tsai, S.-M. et al. Photochemically produced $SO_2$ in the atmosphere of WASP-39b. *Nature* **617**, 483–487 (2023).

22. Leconte, J., Chabrier, G., Baraffe, I. & Levrard, B. Is tidal heating sufficient to explain bloated exoplanets? Consistent calculations accounting for finite initial eccentricity. *Astron. Astrophys.* **516**, A64 (2010).

23. Showman, A. P. & Guillot, T. Atmospheric circulation and tides of “51 Pegasus b-like” planets. *Astron. Astrophys.* **385**, 166–180 (2002).

24. Sainsbury-Martinez, F. et al. Idealised simulations of the deep atmosphere of hot Jupiters. Deep, hot adiabats as a robust solution to the radius inflation problem. *Astron. Astrophys.* **632**, A114 (2019).

25. Sarkis, P., Mordasini, C., Henning, T., Marleau, G. D. & Mollière, P. Evidence of three mechanisms explaining the radius anomaly of hot Jupiters. *Astron. Astrophys.* **645**, A79 (2021).

26. Schneider, A. D., Carone, L., Decin, L., Jørgensen, U. G. & Helling, C. No evidence for radius inflation in hot Jupiters from vertical advection of heat. *Astron. Astrophys.* **666**, L11 (2022).

27. Zhang, K. & Hamilton, D. P. Orbital resonances in the inner neptunian system: II. Resonant history of Proteus, Larissa, Galatea, and Despina. *Icarus* **193**, 267–282 (2008).

28. Millholland, S. Tidally induced radius inflation of sub-Neptunes. *Astrophys. J.* **886**, 72 (2019).

29. Tang, Y. et al. in American Astronomical Society Meeting Abstracts, Vol. 55, id. 151.09 (American Astronomical Society, 2023).

30. Helled, R., Nettelmann, N. & Guillot, T. Uranus and Neptune: origin, evolution and internal structure. *Space Sci. Rev.* **216**, 38 (2020).

31. Anderson, D. R. et al. The discoveries of WASP-91b, WASP-105b and WASP-107b: two warm Jupiters and a planet in the transition region between ice giants and gas giants. *Astron. Astrophys.* **604**, A110 (2017).

32. Welbanks, L. et al. Mass-metallicity trends in transiting exoplanets from atmospheric abundances of H₂O, Na, and K. *Astrophys. J. Lett.* **887**, L20 (2019).

33. Atreya, S. K. et al. in *Saturn in the 21st Century* (eds Baines, K. H., Flasar, F. M., Krupp, N. & Stallard, T.) 5–43 (Cambridge Univ. Press, 2018).

34. Kreidberg, L. et al. A precise water abundance measurement for the hot Jupiter WASP-43b. *Astrophys. J. Lett.* **793**, L27 (2014).

35. Miles, B. E. et al. Observations of disequilibrium CO chemistry in the coldest brown dwarfs. *Astron. J.* **160**, 63 (2020).

36. Menou, K. Turbulent vertical mixing in hot exoplanet atmospheres. *Mon. Not. R. Astron. Soc.* **485**, L98–L103 (2019).

37. Komacek, T. D., Showman, A. P. & Parmentier, V. Vertical tracer mixing in hot Jupiter atmospheres. *Astrophys. J.* **881**, 152 (2019).

38. Mukherjee, S. et al. The Sonora substellar atmosphere models. IV. Elf Owl: atmospheric mixing and chemical disequilibrium with varying metallicity and C/O ratios. *Astrophys. J.* **963**, 73 (2024).

39. Cavalié, T., Lunine, J. & Mousis, O. A subsolar oxygen abundance or a radiative region deep in Jupiter revealed by thermochemical modelling. *Nat. Astron.* **7**, 678–683 (2023).

40. Crossfield, I. J. M. Volatile-to-sulfur ratios can recover a gas giant's accretion history. *Astrophys. J. Lett.* **952**, L18 (2023).

41. Madhusudhan, N. et al. Carbon-bearing molecules in a possible hycean atmosphere. *Astrophys. J. Lett.* **956**, L13 (2023).

42. Chen, J. & Kipping, D. Probabilistic forecasting of the masses and radii of other worlds. *Astrophys. J.* **834**, 17 (2017).

43. Southworth, J. Homogeneous studies of transiting extrasolar planets – IV. Thirty systems with space-based light curves. *Mon. Not. R. Astron. Soc.* **417**, 2166–2196 (2011).

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.

© The Author(s), under exclusive licence to Springer Nature Limited 2024
# Methods

## NIRCam reduction

Both of our new NIRCam observations used the BRIGHT2 readout pattern and the SUBGRISM256 subarray, with the F322W2 observations using seven groups per integration and 1,293 integrations, whereas the F444W observations used 15 groups per integration and 625 integrations. In the subsections below, we will describe the three independent pipelines we used to reduce these NIRCam data.

**Eureka!** Our Eureka! reduction used version 0.9 of the Eureka! pipeline¹⁴, Calibration References Data System (CRDS) version 11.17.0 and context '1093' (for F322W2) or '1097' (for F444W), and jwst package version 1.10.2 (ref. 44). Our reduction methods closely follow those used in previous Eureka! NIRCam spectroscopy analyses¹¹,⁴⁵. The Eureka! Control Files and Eureka! Parameter Files we used are available for download from Zenodo (https://doi.org/10.5281/zenodo.10780448 (ref. 46)) and the important parameters are summarized below.

We ran Eureka!'s Stages 1 and 2 on both the F322W2 and F444W observations using the default jwst processing steps with the exception of increasing the jump rejection threshold in Stage 1 to 6σ from the default of 4σ (to avoid excessive false positives) and turning off the photom step in Stage 2 (as flux-calibrated data are not desired for time-series observations). In Eureka!'s Stage 3, we first carried out two iterations of 5σ clipping along the time axis on the background pixels and we also masked pixels marked as 'DO_NOT_USE' in the data quality array. Next, we corrected for the curvature of the trace by rolling the pixels along the spatial axis in integer pixels (to avoid the need to interpolate or sub-sample the pixels) until the whole trace was approximately centred on the subarray. We then subtracted the background flux from each column in each integration by fitting a linear slope to the column after masking pixels within 13 pixels of the centre of the star's trace and masking any 7σ outliers in the column. The spatial position and point spread function (PSF) width of the star were then recorded for each integration by first summing along the spectral axis and then fitting a Gaussian profile; these parameters were not used in the reduction step but were instead recorded as potential covariates when fitting the light curves in Eureka!'s Stage 5. We then performed optimal spectral extraction⁴⁷ using only the pixels within 5 pixels of the centre of the spectral trace and using a cleaned version of the median integration to compute our spatial profile; to compute this profile, we clipped 5σ outliers along the time axis and then smoothed along the spectral direction using a boxcar filter with a width of 13 pixels. Pixels that differed by more than 10σ compared with the spatial profile were clipped when performing optimal extraction to remove cosmic rays missed during the Stage 1 processing. In Stage 4, we first removed wavelengths with excessively high noise (probably because of unmasked bad pixels); we masked 15 wavelengths in the F322W2 data and one in the F444W data. We then spectrally binned the data into 0.015-µm-wide bins for both NIRCam filters (100 bins spanning 2.45–3.95 µm for F322W2 and 72 bins spanning 3.89–4.97 µm for F444W). Finally, to remove any remaining cosmic-ray effects, we clipped any 4σ outliers in the spectrally binned light curves compared with a cleaned version computed using a boxcar filter with a width of 20 integrations, which was narrow enough to not clip the ingress or egress of the transit.

**Pegasus.** We also reduced and extracted the transmission spectrum of WASP-107b using the Pegasus pipeline (https://github.com/TGBeatty/PegasusProject). We will describe Pegasus in more detail in a forthcoming paper (Beatty et al., in preparation) but we briefly summarize it here. Our Pegasus reduction started from the RATEINTS files from the jwst pipeline v1.10.2, using CRDS version 11.17.0 and context '1093'.

We first performed a background subtraction step for each RATEINTS file by fitting a two-dimensional second-order spline to each integration

using the entire 256 × 2,048 RATEINTS images. We masked out image rows 5 to 75 to not self-subtract light from the WASP-107b system and then we fitted individual background splines to each of the four amplifier regions in the images. We then extrapolated the combined background spline over the masked portions near the star to perform the background subtraction. Visual inspection of the RATEINTS images also showed that, in roughly 5% of the integrations, the reference pixel correction failed for at least one of the amplifier regions, so after the spline fitting and subtraction, we reran the reference-pixel correction using hxrg-ref-pixel (https://github.com/JarronL/hxrg_ref pixels). This seemed to correct the issue.

We then performed a spectrophotometric extraction on the background-subtracted images using optimal extraction techniques⁴⁷. We iteratively constructed a smoothed spatial profile for the F322W2 data using pixel columns from columns 22 to 38 (inclusive) and pixel rows from rows 5 to 1,650 (inclusive), and we iteratively estimated a variance matrix for this same region starting from the pixel uncertainties and read-noise values following established techniques⁴⁷. We then fit a fourth-order polynomial spectral trace to each F322W2 integration and extracted fluxes in each wavelength bin using an extraction aperture with a half-height of 7 pixels centred on the estimated trace in each detector column. In doing so, we accounted for partial pixel effects in both the spectral and the spatial directions.

**tshirt.** We reduced the data using the Time Series Helper & Integration Reduction Tool tshirt¹⁵ (https://github.com/eas342/tshirt), which uses modified steps of the jwst pipeline to calculate rate images and performs a custom covariance-weighted optimal extraction⁴⁸. We followed a very similar procedure as previous NIRCam reductions¹¹,⁴⁵. We begin by running the initial data quality, saturation and superbias subtraction steps with default parameters and reference of the jwst package⁴⁴ version 1.10.2 and CRDS version 11.16.22. For Observation 8 (using the F322W2 filter), CRDS context jwst_1137.pmap, which we reran at a later time. For Observation 9 (using the F444W filter), we context jwst_1093.pmap. Despite the CRDS context numbers, we verified that all reference files were the same between the two filters for Stage 1 processing. We next used a row-by-row odd/even by amplifier (ROEBA) step⁴⁹ to reduce the 1/f/noise by subtracting sky pixels across the spectrum (for the long-wavelength grism) and defocused PSF. We continue the rest of the processing of the calwebb detector1 pipeline with only modification of the jump step's rejection threshold to 6σ. We then use the imaging flat field jwst_nircam_flat_0313.fits and divide the nrcalong images by this for flat-field correction.

For the spectroscopy, the background subtraction was achieved with a column-by-column fit to the regions from 10 to 30 pixels on either side of the source. A robust (sigma-clipped) line was fit to each column of each integration and the line was subtracted from the column. We use covariance-weighted extraction⁴⁸ to sum the pixels in the cross-dispersion direction with an aperture full width of 10 pixels. For the covariance matrix, we assume a constant correlation of 0.08 across pixels. The median measured standard deviation across integration 10 and integration 100 for the F322W2 data with 15-pixel-wide wavelength bins is 1,070 ppm, compared with a theoretical error (from photon and read noise) of 990 ppm. The median measured standard deviation across integration 10 and integration 100 for the F444W data with 15-pixel-wide wavelength bins is 985 ppm, compared with a theoretical error (from photon and read noise) of 971 ppm. We note that F444W has better 1/f subtraction from sky pixels on two amplifiers, whereas F322W2 only has sky pixels available on one amplifier. The noise is higher than theoretical expectations because there are still 1/f correlations in the fast-read direction of the detector.

## HST/WFC3 re-reduction

HST previously observed transits of WASP-107b with WFC3 using the G141 and G102 dispersers. The results of these observations have been
reported previously²,¹⁰. For consistency with the rest of our data analysis, we re-reduced and re-analysed these data.

The G141 data cover a single transit of WASP-107b on UT 6 June 2017 using five orbits of HST time as a part of HST programme GO-14915 (P.I. L. Kreidberg). At the start of each orbit, HST took a single F130N direct image of WASP-107 to establish a wavelength reference and then the observatory switched to spatially scanned spectroscopic observations using the G141 grism. The full details of the observational setup are described in the original report on the results of this observation².

The shorter-wavelength G102 data also cover a single transit of WASP-107b using five orbits of HST time on UT 31 May 2017 as a part of HST programme GO-14916 (P.I.J. Spake). Each orbit of the G102 transit observation began with a F126N direct image to create a wavelength reference and then the G102 grism was used to collect spatially scanned spectroscopy of the WASP-107 system. The complete description of these observations is reported in the original paper describing these data¹⁰.

We reduced both transit observations using a custom data-reduction pipeline that has been described elsewhere⁵⁰ and which we briefly summarize here. For each spatially scanned exposure, we began by extracting 1D spectra for each of the individual up-the-ramp samples, accounting for both the residual background light and the slight angle of the scan pattern on the detector. We then summed these individual spectra to create the 1D spectrum from that exposure. We used the wavelength solutions generated from the direct images taken at the start of each orbit to extract spectroscopic light curves in wavelength bins matching those used in the previous works on these data²,¹⁰. These were nine evenly spaced bins from 0.877 to 1.139 µm for the G102 data and 20 evenly spaced bins from 1.121 to 1.629 µm for the G141 data. We note that, with the revised transit ephemeris for WASP-107b from these observations, we do not see evidence for a starspot crossing in the WFC3 data, as was reported previously², nor do we see starspot crossings in the NIRCam or MIRI observations.

## MIRI re-reduction

We decided to independently re-reduce the MIRI/LRS spectra of ref. 16 to ensure that the reduction was robust to different reduction choices and to ensure a self-consistent reduction of the NIRCam and MIRI data. Our Eureka! re-reduction of the MIRI/LRS observations used version 0.9 of the Eureka! pipeline¹⁴, CRDS version 11.17.0 and context '1097', and jwst package version 1.10.2 (ref. 44). Our reduction method generally follows the Eureka! v1 method described in ref. 51, with some differences and experimental new steps. The Eureka! Control Files and Eureka! Parameter Files we used are available for download from Zenodo (https://doi.org/10.5281/zenodo.10780448 (ref. 46)), and the important parameters are summarized below.

In our Stage 1 processing, we turned on the firstframe and lastframe steps (which respectively discard the first and last frames in each integration) in the jwst pipeline and increased the jump rejection threshold to 7σ from 4σ, as we found that these changes resulted in reduced scatter in the residuals of our light-curve fits. We had also considered different combinations of turning on/off the firstframe, lastframe and Reset Switch Charge Decay (RSCD) correction steps but found that it was best to turn the firstframe and lastframe steps on and leave the RSCD step off for these data. We also investigated using the default CRDS linearity reference file versus the linearity reference file used in ref. 16 and found that there was no clear impact on our final transmission spectra, so we adopted the linearity reference file used in ref. 16.

As pointed out in ref. 51, it seems that narrow wavelength bins result in excessively noisy transmission spectra and underestimated error bars for the MIRI/LRS time-series observation commissioning target L168-9b. Reference 51 proposed that this wavelength-correlated noise could be the result of 390-Hz periodic noise observed in some MIRI subarrays (ref. 52; private communication, Michael Ressler). This noise appears as structured noise with a period of around nine rows within

an individual group (see Extended Data Fig. 1a) and is probably caused by the electronics of MIRI. To investigate the potential impacts of this structured noise on MIRI time-series observations, we developed an experimental 390-Hz noise-removal method for the Eureka! pipeline. This method takes advantage of the knowledge that the noise exhibits a fixed waveform with the phase of the periodic noise remaining fixed within an integration but varying between integrations (only varying slightly and with an apparent cycle of eight integrations); we determined the best-fit waveform using an entire segment of our science observations and then we freely fit the phase shift of the waveform for each integration. We found that there is substantial power in the 1st, 2nd and 4th harmonics (see Extended Data Fig. 1b), with the frequency of the 1st harmonic being exactly 390.625 Hz ($f=1/(256 \times 10 \mu s)$; private communication, Michael Ressler). We found that our new algorithm substantially reduced the periodic noise visible in a Lomb–Scargle periodogram⁵³,⁵⁴ at the group level. While performing this 390-Hz noise-removal step, we also performed group-level background subtraction (GLBS) for each row in each group (the dispersion direction of MIRI is along a column, so rows run in the spatial direction) using the mean value of columns 11–30 and 44–63 (avoiding pixels contaminated by the star); this further removed structured noise from the Lomb–Scargle periodogram. To understand the impact on our final transmission spectra, we also performed another reduction without the GLBS step and another without either the 390-Hz noise removal or the GLBS step.

In Eureka!¹'s Stage 2, we turned off the photom step as flux-calibrated spectra are not desired for time-series observations. In Stage 3, we performed a double-iteration 5σ clipping along the time axis for background pixels, a single iteration of 5σ clipping along the spatial axis for each wavelength and integration-level background subtraction using pixel columns 11–61 (excluding pixels within 10 pixels of the spectral trace) for each wavelength. Similar to the NIRCam spectra, we then performed optimal spectral extraction⁴⁷ using the pixels within 4 pixels of the spectral trace using a cleaned median integration to compute our spatial profile (clipping 5σ outliers along the time axis, smoothing with a 7-pixel-wide boxcar filter and clipping pixels that differed by more than 10σ with respect to the spatial profile). In Stage 4, we first removed three wavelengths that exhibited excessive noise in their light curves. We then binned the spectra using a similar binning scheme as that used in ref. 16 (47 spectral bins spanning approximately 5–12 µm with a constant width of 0.15025 µm); notably, however, we did not use wavelengths below 5 µm, which we expect to be strongly contaminated by light from around 3 µm (ref. 55). We do extend our reduction out to 12 µm as was done in ref. 16 because we do not find evidence for the 'shadowed-region effect' reported by refs. 51,56 in these observations (consistent with the claim that not all observations seem to be affected by this as-yet unexplained phenomenon). Finally, as with NIRCam, we clip 4σ outliers in the spectrally binned light curves compared with a cleaned version computed using a boxcar filter with a width of 20 integrations.

## Derivation of orbital parameters

All of the light-curve fits to the NIRCam, HST and MIRI data described above fix the orbital parameters of WASP-107b to those derived by Murphy et al. (in preparation). We summarize this derivation here.

To precisely determine the orbit of WASP-107b, we simultaneously fit observations of three transits of WASP-107b from the Transiting Exoplanet Survey Satellite (TESS), one from the Goodman High Throughput Spectrograph⁵⁷ (imaging mode, SDSS i-band) on the Southern Astrophysical Research Telescope (SOAR), one from JWST/NIRCam F210M (observed simultaneously with our NIRCam F322W2 data), one from Spitzer/IRAC with IRAC's channel 2 (publicly available from programme 13052, P.I. M. Werner) and one from JWST/MIRI/LRS (broadband, observation described above), as well as radial-velocity observations from CORALIE³¹ and Keck/HIRES⁵,⁵⁸. Together, these transit observations examine a wide range of wavelengths (approximately 0.6–12 µm) and
cover a baseline spanning roughly 6 years, whereas the radial-velocity
observations sample the entire orbit of WASP-107b. We used the default
TESS pipeline reduction, reduced the SOAR/Goodman data using Astro-
ImageJ⁵⁹, reduced the NIRCam/F210M data using tshirt (as described
above), reduced the Spitzer/IRAC data using a custom photometry
pipeline (described in ref. 60) and reduced the MIRI/LRS data using
Eureka! (as described above). For the radial-velocity data, we use the
measurements as tabulated in ref. 5.

We model each of the transits using a batman transit model<sup>61</sup>. For all transit observations except for SOAR/Goodman, we also fit for a background, visit-long linear trend in flux versus time. The SOAR/Goodman data exhibited substantial nonlinear background trends owing to telluric variations during the observation, which we corrected for using a Gaussian process model. We fit all the data described above simultaneously using Markov chain Monte Carlo (MCMC) sampling with emcee (ref. 62), sampling the time of conjunction, orbital period, inclination, semimajor axis, eccentricity, argument of periastron, planet-to-star radius ratios in each bandpass, quadratic limb-darkening coefficients in each bandpass, radial-velocity semiamplitude and system velocity and the slope and intercept of the linear background trend of each visit. We ran this sampling for 10,500 steps, which was sufficient for each parameter to converge (&gt;25× the average autocorrelation times, plus we visually inspected each MCMC walker time series). The best-fit parameter values are given by Murphy et al. (in preparation) and were P = 5.72148722 ± 3 × 10<sup>-7</sup> days, t<sub>0</sub> = 2,459,958.747244 ± 8.2 × 10<sup>-6</sup> BJD<sub>TDB</sub>, i = 89.57 ± 0.03°, a/R<sub>*</sub> = 18.05 ± 0.1, e = 0.05 ± 0.01 and ω = -2.3 ± 6.1°.

Light-curve fitting

**Eureka!** In Eureka!'s Stage 5, we fitted astrophysical and systematic noise models to our spectroscopic light curves. For both NIRCam filters, our systematic noise model consisted of a quadratic trend in time and a linear decorrelation as a function of the spatial position and PSF width measured during Stage 3. For MIRI/LRS, we used a linear trend in time, an exponential ramp in time and a linear decorrelation as a function of the spatial position and PSF width. For all observations, our astrophysical model consisted of a starry<sup>63</sup> transit model with quadratic limb-darkening parameters fixed to an ATLAS model limb-darkening spectra<sup>64</sup> as computed by the Exoplanet Characterization Toolkit (ExoCTK)<sup>65</sup>. We fixed all of our orbital parameters to those of Murphy et al. (in preparation; see the 'Derivation of orbital parameters' section above) and adopted a minimally informative prior on the planet-to-star radius ratio.

Following the recommendations in ref. 51, we removed the first 800 integrations of the MIRI/LRS observations to remove the worst part of the initial ramp in the light curve. We also removed integrations 475–490 from the NIRCam F444W observations that exhibited a brief spike in noise. We removed no further integrations from the NIRCam F322W2 observations. Finally, we also fitted a white-noise multiplier to account for excess white noise. For all our observations, we used PyMC3’s No U-Turn Sampler<sup>66</sup> using two independent chains, each taking 4,000 tuning draws and then 3,000 posterior samples with a target acceptance rate of 0.85. We then validated that the chains had converged by ensuring that the Gelman–Rubin statistic<sup>67</sup> was at or below 1.1. The 16th, 50th and 84th percentiles of the posterior samples were then used to estimate the best-fit values and uncertainties for each fitted parameter.

For the NIRCam F322W2 observations, we find that the white noise is,
on average, about 15% above our estimate of the stellar photon-limited
noise floor in each channel, whereas F444W exhibited white noise only
about 7% above the estimated photon limit on average. This excess noise
is probably the result of insufficiently corrected 1/f noise<sup>48</sup>, which gets
partially converted to white noise in the time series when spectroscopi-
cally binning the light curves. For MIRI/LRS, we used the estimated gain
of 3.1 electrons/DN from refs. 51,56, as the current value of 5.5 electrons/
DN used in the CRDS reference files is known to be incorrect. With this

estimated gain, we find that our fitted white-noise level is only about
10% above the estimated stellar photon noise limit near 5 µm but climbs
steadily to approximately 50% above the limit near 10 µm and then
climbs steeply to around 140% above the limit near 12 µm. Part of this
excess long-wavelength white noise is probably the result of the impact
of background noise, as these limits only considered the stellar photon
noise. Although neither of our NIRCam observations show evidence for
residual red noise in their Allan variance plots<sup>68</sup>, several of our MIRI/LRS
channels do show moderate amounts of residual red noise. Following
the β error-inflation method developed in ref. 69, we computed β, the
ratio of the median value of our Allan variance plots at 15–23-min time-
scales (around the approximately 19-min transit ingress/egress dura-
tion of WASP-107b) to the value expected for white noise and multiplied
our transmission spectra uncertainties by this amount. This resulted in
an increase in the uncertainties by about 50% for wavelengths ≤6.5 µm
and 0–30% for longer wavelengths. After this, we arrived at our final
2.45–12.00 µm transmission spectra.

**Pegasus.** We fit the F322W2 and F444W spectroscopic light curves from our Pegasus reduction using a batman transit model⁶¹. We fixed most of the transit-model parameters to the values measured from broadband fits to the JWST data (Murphy et al., in preparation; see the 'Derivation of orbital parameters' section above), leaving only the planet-to-star radius ratio as the only free astrophysical parameter. We also simultaneously fit for a background quadratic trend across each visit, which had two associated slope parameters and a normalization factor. We did not impose a prior on any of these four parameters and we fit each spectroscopic channel individually. For each channel, we used a quadratic limb-darkening law and fixed the limb-darkening coefficients to the values calculated from the ATLAS model⁶⁴ by ExoCTK⁶⁵ in the F322W2 and F444W bandpasses. To generate the limb-darkening coefficients, we used the spectroscopic stellar properties measured in a previous work⁵.

To perform the light-curve fit in each spectroscopic channel, we
performed an initial Nelder-Mead likelihood maximization followed
by MCMC likelihood sampling. We initialized the MCMC chains about
the maximum likelihood point estimated from the Nelder-Mead maxi-
mization. We used the emcee⁶² Python package to perform the MCMC
sampling, using 20 walkers with a 2,000-step burn-in and a 4,000-step
production run. At the end of this process, we checked that the
Gelman-Rubin statistics was below 1.1 for each parameter in each
spectral channel to ensure that the MCMC sampling had converged.

We performed two tests to evaluate the goodness-of-fit of our transit modelling in each spectral channel. First, we confirmed that the per-point flux uncertainties in the light curves were consistent with the standard deviation of the best-fit model residuals. Second, we calculated the Anderson-Darling statistic for the residuals to verify that they followed a Gaussian distribution. Our spectroscopic light-curve fits did not show any evidence of non-Gaussianity in the residuals.

The estimated depth uncertainties were 1.15× the photon-noise expectation in the F322W2 data and 1.06× the photon-noise expectation in the F444W data.

tshirt. We followed a similar procedure as previous NIRCam light-curve
fitting<sup>11,45</sup>. We use the starry<sup>63</sup> code to model the light curves and PyMC3
(ref. 66) to sample the posterior distributions of the fits. We adopted
uninformative priors on the stellar limb darkening with a two-parameter
quadratic law<sup>70</sup>. This allows for the comparison of results with the
ATLAS stellar limb-darkening models. The astrophysical starry model of
the light curve is multiplied by a second-order (quadratic) polynomial
baseline to allow for trends from instrument and astrophysical stellar
variability. We bin the spectra in the time axis to 300 bins for faster
light-curve evaluation with starry. Finally, we include a parameter for
the standard deviation of the data to allow for excess noise beyond
the theoretical photon and read noise. We sampled the posterior with
the No U-Turn Sampler<sup>66</sup>. We used 3,000 tuning steps, 3,000 sampling steps and two chains with a target acceptance rate of 90%.

**Comparing NIRCam reductions.** A motivation for conducting three different reductions of the NIRCam data was to check the robustness of our final transmission spectrum against different choices and assumptions made during the image calibration, spectral extraction and light-curve fitting stages of our analyses. Generally, we find that all three reductions give the same general transmission spectrum with approximately the same depth uncertainties (Fig. 2 and Extended Data Fig. 2a). The mean depth uncertainty across the entire NIRCam wavelength range is approximately 90 ppm at our constant wavelength $\Delta\lambda = 0.015$ µm binning, which is larger than the mean point-to-point difference between the Eureka! and Pegasus spectra (58 ppm) and the Eureka! and tshirt spectra (72 ppm). A $\chi^2$ comparison of the three spectra gives $\chi^2/\text{dof} = 0.74$ for Eureka! versus Pegasus and $\chi^2/\text{dof} = 0.64$ for Eureka! versus tshirt. The larger difference between Eureka! and tshirt is mainly because the F322W2 spectrum of tshirt consistently falls below that of Eureka! and Pegasus at wavelengths longer than about 2.893 µm; this difference may be caused by the ROEBA method of tshirt because the deviation occurs exactly at an amplifier boundary (see Extended Data Fig. 2a). In both pipeline comparisons, the small reduced chi-squared values imply a $<0.1\sigma$ likelihood that the spectra are drawn from two different underlying distributions for our 172 degrees of freedom. Because the choice of reduction methods did not substantially affect our final NIRCam spectra, we chose the Eureka! reduction as our fiducial NIRCam spectra for all subsequent atmospheric modelling work.

To further investigate the sensitivity of our results to different limb-darkening assumptions, we also ran four separate fits with our fiducial Eureka! pipeline, namely, fixing our limb-darkening parameters to: (1) ATLAS-predicted quadratic law coefficients computed using ExoCTK (refs. 64,65); (2) STAGGER-predicted quadratic law coefficients computed using ExoTiC-LD (refs. 71,72); (3) STAGGER-predicted four-parameter law coefficients computed using ExoTiC-LD (refs. 71,72); and (4) freely fitting quadratic law coefficients using the reparameterization of ref. 70 for efficient and physical sampling. As shown in Extended Data Fig. 2b, none of these spectra markedly differ from each other. In particular, we find that comparing the transmission spectrum of our fiducial analysis that used the ATLAS quadratic limb-darkening coefficients to fits with limb-darkening choices only results in a reduced chi-squared value of 0.02 for STAGGER's quadratic model, 0.09 for STAGGER's four-parameter model and 0.28 for the fit with freely fit reparameterized quadratic limb-darkening coefficients. In all three cases, this means that the differences between the spectra are much smaller than the uncertainties on the measured depths. We also further disfavour the freely fit limb-darkening coefficients, as these fitted coefficients exhibited substantially more point-to-point scatter than would be expected for a physical limb-darkening spectrum.

**HST/WFC3.** We fit each of the 29 G102 and G141 spectroscopic transit light curves using a batman transit model<sup>61</sup> and standard detrending techniques for WFC3 time-series data<sup>50</sup>. As with the JWST spectroscopic data, for these fits, we fixed the transit centre time, orbital period, orbital inclination, scaled semimajor axis, orbital eccentricity and the argument of periastron to previously measured values (Murphy et al., in preparation; see the ‘Derivation of orbital parameters’ section above). We also fixed the quadratic limb-darkening coefficients in each spectroscopic channel to the values estimated from the ATLAS limb-darkening model<sup>64</sup> by the ExoCTK<sup>65</sup>. This left the planet-to-star radius ratio as the only free astrophysical parameter. We fit each spectroscopic channel individually.

All of the G102 and the G141 light curves showed the usual ‘fishing-hook’ systematic trend within each orbit and the background linear slope typical for WFC3 time-series observations. We modelled this in

each spectral channel as a part of the transit fitting process using an exponential ramp within each individual orbit, and a linear background trend across each visit, of the form

$$F_{\text{detrend}} = (mt_{\text{V}} + n) \left( 1 - Ae^{\frac{t_{\text{O}}}{\tau}} \right). \qquad (1)$$

Here $t_O$ is the time since the start of observations of each orbit, $t_V$ is the time since the start of the visit, A and $\tau$ are the fitting coefficients for the exponential ramp and m and n are the fitting coefficients for the linear trend. As well as fitting the systematics in this way, we also did not use the data from the first orbit within each visit nor the first exposure from each orbit when fitting the data.

To fit the transit light curve in each spectroscopic channel, we first performed a Nelder-Mead likelihood maximization to find an initial best fit. We then used the emcee (ref. 62) Python package to perform a MCMC exploration of the surrounding likelihood space to improve this best-fit estimate. In each channel, we used 20 MCMC walkers, each with a 1,000-step burn-in followed by a 4,000-step production run. At the end of this, we judged the MCMC to have converged by verifying that the Gelman-Rubin statistic for each parameter was less than 1.1. We also checked to ensure that the median per-point flux uncertainty matched the standard deviation of the residuals to the best-fit light-curve model in each channel. We also performed an Anderson-Darling test on the residuals of each channel to check for non-Gaussianity. The light-curve residuals in each channel seem well behaved.

We note that we do not see evidence for a starspot crossing in the third orbit of the G141 transit data, as was suggested previously<sup>72</sup>. Instead, the revised transit ephemeris and orbital properties we measure using the JWST data for WASP-107b serve to shift the start of transit egress slightly earlier at the time of the G141 observations. This accounts for the putative starspot feature in the residuals of the earlier analyses. As a result, we do use the data from the third G141 orbit in our fitting. This slightly improves our measured depth uncertainties compared with those of previous results.

**MIRI.** Unfortunately, although our 390-Hz noise-removal and GLBS steps remove structured noise from the Lomb-Scargle periodogram<sup>53,54</sup> of the MIRI/LRS group-level data, we find that these steps ultimately have a fairly small impact on the final transmission spectra of our science target WASP-107b, as well as the MIRI/LRS time-series observation commissioning target L168-9b. For example, Extended Data Fig. 3 shows that the MIRI/LRS observations of L168-9b continue to show excess noise after the 390-Hz noise removal step has been applied. It is possible that the 390-Hz noise-removal step does eliminate excess noise for 1-µm bin sizes, but at that coarse a resolution, there are only seven bins and there is substantial uncertainty in the estimated standard deviation from the small sample size. Also, with finer bin-size sampling than was used in ref. 51, we are better able to understand how the excess noise varies with varying bin size. In particular, whereas we find that the excess noise decreases with increasing bin size (similarly to what was reported in ref. 51), there also seems to be a spike in the excess noise around a bin width of 0.25 µm; this corresponds to an average spectral bin width of 9 pixels (although the bin width varies with wavelength), which is also the approximate period of the 390-Hz noise. Surprisingly, although this noise spike seems to share an approximate period with the known 390-Hz noise, the 390-Hz noise-removal step does not seem to have removed the excess noise in the final transmission spectrum of L168-9b.

We also see very limited impact from both the 390-Hz noise-removal and GLBS steps in the WASP-107b MIRI/LRS science observations (see Extended Data Fig. 4). Most probably, this is because there is fairly minimal curvature in the waveform over a single row (which only spans about 7% of the first harmonic, which has the largest amplitude), so the per-row background-subtraction step performed on the
integration-level data in Stage 3 adequately subtracts this structured noise. It is possible that this 390-Hz noise-removal procedure could still be useful to future works that seek to use separate background calibration observations taken before and/or after the science observations instead of computing the background on the science observations themselves, but—at present—we do not recommend the 390-Hz noise-removal procedure be applied to science observations alone owing to the minimal impact and high computational cost. That said, it does not seem that our 390-Hz noise-removal procedure negatively affects our results and we ultimately decide to use it in our fiducial reduction. The impacts of the GLBS step are noticeable but also minor, and it is as-yet unclear whether the noise introduced by estimating the group-level background outweighs the structured noise removed by this step.

Comparing our fiducial spectrum (with the 390-Hz noise-removal and GLBS steps) with the European Consortium's three reductions presented in ref. 16 (see Extended Data Fig. 4), which used the CASCADe and TEATRO pipelines, as well as the Eureka! pipeline, it is again clear that the overall shape of the spectrum is robust to different analysis pipelines and different reduction choices. As pointed out in ref. 16, the transit depth differences among the European Consortium's spectra show an increasing trend with wavelength, with the CASCADe pipeline giving larger depths at longer wavelengths, which they attribute to different systematics models when fitting the light curves. As shown in Extended Data Fig. 4, we find fairly similar conclusions with our fiducial spectrum, which generally agrees well with the CASCADe spectrum. Looking more closely, it seems that our fiducial spectrum has a slightly smaller transit depth around 7.5 µm and slightly larger transit depths around 7 µm and from around 8–10 µm compared with the three reductions in ref. 16. This results in our fiducial spectrum having a slightly smaller bump near the $v_3$ SO₂ feature and a slightly larger bump near the $v_1$ SO₂ feature.

## Atmospheric retrievals

Detailed interpretation of atmospheric spectra requires comparisons with atmospheric radiative transfer models by means of statistical algorithms such as Bayesian inference methods. These data-model inference techniques are commonly known as 'atmospheric retrievals' and enable constraints on the atmospheric properties of an exoplanet, such as chemical composition and vertical temperature structure from an observed spectrum. The inference framework computes transmission spectra, generally on the order of millions or tens of millions per inference, from a parametric atmospheric model (that is, a 'forward model') that solves line-by-line radiative transfer under hydrostatic equilibrium. These atmospheric models generally assume a plane-parallel atmosphere and include parameters that determine the chemical abundances for different chemical gases, the vertical pressure–temperature structure of the planet and the presence of clouds/hazes as further sources of opacity in the atmosphere. Below we present the two atmospheric modellings used in this study, spanning a wide range in physico-chemical assumptions.

**1D-RCPE grid retrieval.** To produce physically self-consistent solutions, we fit the transmission spectrum with a suite of models that impose 1D-RCPE. The coupling between the radiative–convective equilibrium solver (ScCHIMERA<sup>11,73,74</sup>) and the kinetics code (VULCAN<sup>21,75</sup> using the H–O–C–N–S reaction list from ref. 21), as well as further details, has been previously described in ref. 11. The coupled model requires as inputs the incident stellar flux spectrum ($T_{\text{eff}}$ = 4,425 K, log(g) = 4.63 c.g.s. (ref. 76)) at the top of the planetary atmosphere (including the ultraviolet for photochemistry, using the HD 85512 MUSCLES spectrum<sup>77–79</sup> as a proxy, which matches the effective temperature and gravity of WASP-107 well within 1%), which can be scaled through an effective irradiation temperature, an internal temperature to set the deep adiabat, the atmospheric elemental abundances (through metallicity and a C/O as

described in ref. 11), an eddy diffusivity and the bulk planet properties, radius and mass. The outputs are the converged 1D temperature and gas volume mixing ratio pressure profiles that can then be 'post-processed' to produce a transmission spectrum. Extended Data Fig. 5 highlights select parameter slices through the grid and the subsequent impacts on the spectra.

As this process is computationally expensive owing to the large numbers of converged models required for the Bayesian inference, we generate the grid in stages. First, we generated a course grid in metallicity ($0 \le [M/H] \le 2$ in steps of 0.5), C/O ($0.1 \le C/O \le 0.7$ in steps of 0.2), internal temperature ($100 \le T_{\text{int}} \le 500$ in steps of 100 K) and a vertically constant eddy mixing ($7 \le \log_{10}K_{zz} \le 9$ in steps of 1) but at a fixed effective irradiation temperature (at the planetary equilibrium temperature, 738 K). We then perform Bayesian inference on the observed spectrum with this course grid as described in ref. 11 (with more details specific to WASP-107b described below). In the fitting process, we also included a temperature profile offset free parameter (a simple additive shift to the whole profile) to account for possible variations in temperature between the limb and the planetary dayside, which—through the scale height—can influence the feature sizes.

The results of the course-grid exercise effectively 'narrowed' the plausible parameter space down, at which point we generated a more finely sampled grid (below which the grid spacing had no effect). This 'fine' grid was generated over [M/H] ($0.5 \le [M/H] \le 1.875$ in steps of 0.125), C/O ($0.05 \le C/O \le 0.6$ in steps of 0.05), internal temperature ($200 \le T_{\text{int}} \le 550$ in steps of 50 K) and eddy mixing ($7 \le \log_{10}K_{zz} \le 10$ in steps of 0.5), resulting in 6,048 converged 1D-RCPE models. On the basis of the temperature offset parameter from the course grid fit, this fine grid was generated at a cooler irradiation temperature (560 K) to more appropriately adjust for changes in the chemistry with temperature. This finer grid (at the new fixed irradiation temperature) was then refit to the data using the same inference procedure and is what we use to inform our primary conclusions. We also tested the effect of temperature by fitting, again, for a temperature profile offset parameter and found that it was consistent with zero at nearly $10$ ($-18 \pm 14$ K), suggesting that our choice of irradiation temperature for the fine grid resulted in a temperature profile that was consistent with the observed spectrum.

In the Bayesian inference procedure, as described in ref. 11, we also fit for, on-the-fly, the planetary reference pressure and reference radius (ref. 80, which shifts and stretches the planetary spectrum and also affects the zero point for the planetary gravity with height), an offset for the MIRI spectrum relative to the NIRCam F444W spectrum, the temperature profile offset parameter (described above), an ammonia abundance enhancement factor, a power law haze + grey cloud and a parametric cloud profile to describe unknown cloud resonance features over the MIRI wavelengths (described more below). Including these and the 1D-RCPE grid parameters results in a total of 17 free parameters. The transmission spectra are generated with $R = 100,000$ absorption cross-sections (same line lists cited in ref. 11 but now also including SO₂ (ref. 81)) and fit to the data with the PyMultiNest routine, using a total of 500 live points and the default parameters. Gas and cloud detection significances are computed as described in refs. 32,82. Generally, the detection of different gases in the 1D-RCPE solution corresponds to a preference for the opacity of the gas in question as a substantial contribution to the spectrum, because only the opacity is removed from the inference and not from the chemistry, as described in ref. 11. The retrieved parameters and priors are shown in Extended Data Table 1.

To illustrate the sensitivity of our results to the key physical processes that enabled us to arrive at the hot-interior solution, we explored perturbations to the atmospheric structure and resultant spectra, shown in Extended Data Fig. 5. We perturb the nominal 1D-RCPE solution (given in Extended Data Table 1) by: (1) turning off mixing and photochemistry (equilibrium versus disequilibrium, top row of Extended Data Fig. 5); (2) changing the internal/effective temperature (middle row of Extended
Data Fig. 5); and (3) changing the strength of the eddy diffusion (bottom row of Extended Data Fig. 5). Mixing from a hot deep layer (arising from a high $T_{int}$) results in an approximately three-orders-of-magnitude depletion in the methane abundance compared with equilibrium. This has a substantial influence on the spectrum and readily explains the 'muted' methane features. Equilibrium abundances are clearly ruled out by the free-retrieval constraints ($-6.0 < \log_{10} \text{CH}_4 < -5.6$; details below). It is also apparent that $T_{int}$ values lower than about 350 K and a $\log_{10} \text{K}_{zz}$ lower than about 8 struggle to result in the necessary reduction in the methane abundance, even with mixing.

**Cloud parameterization.** There are two leading models for modelling clouds and hazes in exoplanetary atmospheres: microphysical models and parametric models. The former requires understanding of the nucleation, condensation, evaporation, coagulation and transport of clouds<sup>83</sup>. The latter aims to capture the spectroscopic signature of these condensates and aerosols without any assumptions of the physical and chemical processes that lead to their formation and destruction, attempting to account for their presence to unbias any inferences of other properties of interest (for example, volume mixing ratios of different gases or atmospheric metallicity; refs. 18,84–86). Some parametric models have attempted to capture the functional dependence of cloud extinction on wavelength using analytic models<sup>87,88</sup> or Mie theory<sup>89–92</sup>.

Of relevance to this work, one of the approaches in ref. 16 to fit the MIRI + HST/WFC3 G141 spectra of WASP-107b was to use the eddysed microphysical cloud parameterization<sup>93,94</sup> for several silicate condensates. Briefly, this cloud framework assumes a balance between upwards turbulent mixing (parameterized with a vertical eddy diffusivity, $K_{zz}$) and sedimentation (parameterized with $f_{sed}$) of droplets. This balance governs the particle-size distribution and abundance with height in the atmosphere (larger droplets deeper towards the cloud base and smaller droplets at higher altitudes). The overall abundance of each cloud species is set by the condensate abundance at the base of the cloud. The droplet abundances/profiles along with condensate/droplet optical properties (derived from Mie theory) are then used to compute the total cloud extinction. Specifically, they fit for a cloud-base pressure, condensate abundance, $f_{sed}$ and $K_{zz}$ assuming no self-consistency in the location and abundance of the condensates with the chemistry or temperature profile. They found that $\text{SiO}_2$ best explained the MIRI spectral shape. However, the retrieved location of the cloud base was at relatively low pressures (on the order of mbar)—at much cooler temperatures and lower pressures than otherwise would be expected from equilibrium condensation chemistry (Extended Data Fig. 6a)—and a highly compact cloud (large droplets).

We initially follow the same eddysed cloud-parameterization approach described above (assuming Mie theory with the indices of refraction from ref. 95). Although we are able to find satisfactory fits to the MIRI + HST observations (as in ref. 16), when applied to the entire broadband spectrum with our NIRCam observations, we were unable to find satisfactory fits.

We also tested the plausibility of silicate clouds (MgSiO<sub>3</sub>, SiO<sub>2</sub> and Mg<sub>2</sub>SiO<sub>4</sub>) and salt/sulfide clouds (Na<sub>2</sub>S and KCl) in explaining the broadband spectrum within the same eddysed<sup>93,96</sup> framework. We use the 1D-RCPE atmospheric structure from an initial representative fit (temperature-pressure profile and chemistry using the nominal composition [M/H] = 1.0, C/O = 0.3, eddy diffusion coefficient $\log K_{zz} = 8.5$ and $T_{int} = 500$ K) to self-consistently determine which clouds and where they form along the temperature profile. We varied $f_{sed}$ between 0.1 (vertically extended clouds) and 1.0 (compact clouds). Extended Data Fig. 6 summarizes the condensate mixing ratios (panel a) and their impact on the spectrum (panel b) for $f_{sed} = 0.2$. Only the lower $f_{sed}$ values (<0.5, more extended clouds) are able to produce a notable impact on the spectrum, as higher values result in clouds that are too compact with much of the cloud opacity at too deep pressures. However, the lower $f_{sed}$ values result in smaller droplets, which produce a more narrow

resonance feature than what is needed to fit the MIRI observations. It is the combination of a poor fit to the full broadband spectrum along with the implausibility of the silicate cloud bases to exist at the altitudes revealed by the observations that instead prompts us to take a phenomenological approach to the cloud modelling.

We introduce a new parametric treatment to capture the spectroscopic signatures of cloud particulate resonance feature. This approach is motivated by previous efforts to model optical slopes resulting from hazes as a scaling to Rayleigh scattering<sup>18</sup>. We model these concave spectroscopic signatures (see Fig. 2 of ref. 16) of cloud condensates using a Gaussian function (that is, the composition of an exponential function and a concave quadratic function). The total extinction coefficient is given by a skewed normal distribution of the form

$$ \kappa_{\text{cond}}(\lambda) = 2\kappa_{\text{opac}}\phi\left(\frac{\lambda - \lambda_0}{\omega}\right)\phi\left(\frac{\xi - \lambda - \lambda_0}{\omega}\right) \quad (2) $$

in which $\kappa_{\text{opac}} = K_{\text{cond}}n_{\text{tot}}$ is the extinction coefficient for an opacity free parameter $K_{\text{cond}}$ and $n_{\text{tot}}$ is the total number density of the atmosphere, $\phi$ is the standard normal probability density function and $\Phi$ is the standard normal cumulative density function, $\lambda_0$ is the wavelength at which the Gaussian is centred (that is, the mean of the Gaussian), $\omega$ is the scale parameter (that is, the standard deviation of the Gaussian) and $\xi$ is a shape parameter that controls the skewness of the distribution so that, when $\xi = 0$, the standard normal distribution is recovered.

In both the 1D-RCPE retrievals and free retrievals, we incorporate this parametric treatment alongside scalings to the Rayleigh scattering in the optical<sup>18</sup> and optically thick grey clouds, as explained in ref. 11, in a linear combination with a cloud-free model to account for the presence of cloud/haze inhomogeneities<sup>84</sup>. Future studies may investigate the best practices for the parameterization of cloud particulate resonance features.

**Free retrieval.** We further explore the atmospheric properties inferred from the spectrum of WASP-107b using more flexible and agnostic forward models in a Bayesian inference procedure. These, known as 'free retrievals', are methods to retrieve the chemical abundances of different gases, the vertical temperature structure and the cloud/haze properties on the planetary atmosphere. Compared with the 1D-RCPE grid retrieval explained above, these methods do not assume any physico-chemical equilibrium conditions and instead aim to capture the atmospheric conditions directly through a series of parameters for the chemistry and physical conditions of the atmosphere without expectations of physical consistency. These more flexible approaches provide the opportunity to capture conditions that otherwise would be prohibited by self-consistent models, such as combinations of gases not considered under chemical equilibrium. Nonetheless, caution must be exercised in the interpretation of these free retrievals as model assumptions may contribute to biased and unphysical atmospheric estimates (see ref. 85 for a discussion). As with the 1D-RCPE grid retrieval above, we simultaneously retrieve on the HST/WFC3 G102 and G141, JWST/NIRCam F322W2 and F444W and JWST/MIRI observations.

We use two independent free-retrieval frameworks in our analysis: Aurora<sup>91</sup> and CHIMERA<sup>2,84,97</sup>. We select the latter as our fiducial free-retrieval comparison because the radiative transfer code, sources of opacity and model resolution is the same as in the 1D-RCPE grid retrieval above, therefore enabling a more direct comparison. The former, Aurora, is used to consider the impact of different sources of opacity than those used by CHIMERA and consider the impact of line-by-line cross-section at higher sampling resolution. We find that our detection significances and conclusions are largely independent of the framework used, with the details of each analysis described below.

CHIMERA solves radiative transfer for a parallel-plane atmosphere, under hydrostatic equilibrium, for transmission geometry. The atmospheric model considers a one-dimensional model atmosphere spanning
from about $10^{-9}$ bar to about 100 bar, divided into 100 layers uniformly spaced in logarithmic pressure space. The vertical temperature structure is parameterized following the prescription in ref. 98. The model simultaneously retrieves the reference pressure and reference radius for the assumed planetary gravity of $log_{10}g = 2.45$ c.g.s.

The atmospheric models assume uniform mixing ratios for H₂O, CH₄, NH₃, CO, CO₂, SO₂ and H₂S using independent free parameters for the volume mixing ratio of each gas. Our choice of chemical species is motivated by those expected in exoplanetary atmospheres at these warm temperatures⁹,¹⁷. As part of our initial analysis, we considered the presence of PH₃, as this species may be expected for some substellar objects at these temperatures. Nonetheless, the observations did not place meaningful constraints on the abundance of PH₃ and was not considered in our fiducial run in an effort to limit the number of free parameters in our analysis. The presence of clouds and hazes is considered by using the one-sector parameterization in ref. 91 introducing the combined spectroscopic effect of a optically thick cloud deck with a cloud opacity $\kappa_{cloud}$ and hazes following an enhancement to Rayleigh scattering¹⁸ in a linear combination with a cloud-free atmosphere⁸⁴. Furthermore, we incorporate the effect of unknown cloud resonances (for example, wavelength-dependent condensates) as explained above.

The Bayesian inference is performed using nested sampling⁹⁹,¹⁰⁰, through PyMultiNest¹⁰⁰⁻¹⁰², using 500 live points in the sampling. Each forward model in the sampling was calculated using line-by-line opacity sampling at a spectral resolution of 100,000 and then binned to the resolution of the observations. The line lists considered remain the same as in the 1D-RCPE analysis described above. In total, the sampling is performed over 24 parameters: seven molecular gases, six for the pressure-temperature structure of the planet, one for the reference pressure, one for the reference radius, eight for the presence of inhomogeneous clouds and hazes and one for an offset between the JWST/MIRI observations and the JWST/NIRCam F444W observations.

Our second free retrieval with Aurora⁹¹ follows largely the same model setup as that described above with CHIMERA. A detailed description of the transmission modelling approach of Aurora is available in ref. 91. The main differences relative to the retrieval with CHIMERA are: (1) different opacity sources for H₂O, CO₂, NH₃, CH₄ and H₂S obtained from refs. 81,103–105; (2) the use of a free parameter $P_{cloud}$ to account for optically thick clouds instead of the $\kappa_{cloud}$ parameter; and (3) computing the forward models using line-by-line cross-section sampling at a spectral resolution of 20,000 instead of 100,000. The retrieved atmospheric properties are generally consistent within the CHIMERA retrieval at 68% confidence and with predictions from the 1D RCPE grid retrieval. The results from this retrieval are also included in Extended Data Table 2. Supplementary Fig. 2 shows the equivalent to Fig. 4 with the Aurora constraints and highlights that, although the inferences are consistent, the use of lower spectral resolution results in wider constraints in this case.

Supplementary Fig. 1 shows the retrieved transmission spectrum and the contributions of the detected gases in the atmosphere of WASP-107b as inferred with Aurora. In both cases, the retrieved volume mixing ratios from the free retrievals are consistent with the gas profiles from the 1D-RCPE inferences, as shown by the posterior distributions from the free retrieval and gas samples from the 1D-RCPE model in Fig. 4 and Supplementary Fig. 2. Furthermore, the retrieved vertical pressure-temperature structure of the planet is consistent within the fully parametric free retrieval and the self-consistent 1D-RCPE model as shown in Extended Data Fig. 7. The free retrieval with Aurora finds an offset between the JWST/MIRI and JWST/NIRCam F444W observations of $282^{+49}_{-41}$ ppm. A complete summary of the retrieved parameters and their priors is shown in Extended Data Table 2.

Both free retrievals confirm the detections of H₂O, CH₄, NH₃, CO, CO₂ and SO₂ with detection significances of 5σ or greater. The reported detection significances in the main text result from the model

comparisons with CHIMERA, whereas the detections with Aurora remain equally substantial at 5σ or greater. Similarly, we find that the fiducial model with optically thick clouds, Rayleigh-enhancement hazes and wavelength-dependent condensates is preferred over a cloud-free atmosphere at 26σ. Including wavelength-dependent condensates is preferred at 13σ compared with just including optically thick clouds and Rayleigh-enhancement hazes. We note that the term 'detection significance' refers to a model preference between a reference model considering all gases and nested models for which each individual species is removed in turn⁸²,⁹¹. As such, any quoted detection significance is dependent on the choice of models and choice of model priors. Moreover, the conversion between differences in Bayesian evidence and 'sigma detections' can result in extremely large values (for example, >10σ) that are difficult to interpret (see ref. 106 for a discussion). For WASP-107b, the agreement in retrieved molecular abundances and confirmation of the strong evidence for the detection of these species regardless of model assumptions (for example, free abundances versus radiative-convective-photochemical equilibrium) confirms the robustness of the interpretation of the planet's spectrum.

**Robustness of ammonia detection.** To provide scrutiny beyond using Bayesian evidence model selection for the NH₃ detection, we perform leave-one-out cross-validation (LOO-CV (refs. 107–109)) on the models with and without NH₃ from the 1D-RCPE grid retrievals. In LOO-CV, a model is trained on the dataset with one data point left out and then scored according to how it can predict the left-out data point (that is, by calculating the expected log predicted density of the left-out data point, elpd). This procedure is repeated for all data points allowing the out-of-sample predictive performance of a model to be estimated. The LOO-CV scores for each data point can then be compared between two models, highlighting where one model outperforms the other¹⁰⁶,¹¹⁰,¹¹¹ and, consequently, in this case, which wavelengths and instruments drive the NH₃ detection. Extended Data Fig. 8 shows the difference in elpd (Δelpd) per data point in the spectrum of WASP-107b for the reference 1D-RCPE model and the model without NH₃ absorption.

The Bayesian model comparison finds a 6σ detection (that is, model preference) for the model including NH₃ absorption. By performing a LOO-CV analysis of these models, we can determine which data points drive this model preference. We find that the inclusion of NH₃ improved the model performance at all wavelengths (for example, there are points with positive Δelpd at all wavelengths), with a localized region of high-performance data points near 3 µm. Indeed, the points with the highest preference for the model with NH₃, those with Δelpd > 1, are all between 2.9 and 3.2 µm and part of the NIRCam F322W2 observations. After NIRCAM F322W2, the next instrument with positive Δelpd scores is HST-WFC3 G141 (9th highest scoring point), followed by HST-WFC3 G102 (14th highest scoring point) and NIRCAM F444W (17th highest scoring point). The highest scoring point of MIRI is at 10.17 µm and corresponds to the 25th highest scoring point with Δelpd ≈ 0.4.

We proceed to total the difference in LOO-CV scores (that is, $sum(\Deltaelpd)$) between the models with and without NH₃ over regions in which the NH₃ cross-section is dominant (that is, for which NH₃ contributes >50% of the total cross-section). We find that the density of the increased predictive performance (that is, $sum(\Deltaelpd)/\text{number of points}$) is higher, and more than double the value, in regions in which NH₃ is dominant than in regions in which NH₃ is not the dominant cross-section. This confirms that the detection of NH₃ is greatly improved by the data in which substantial absorption by the gas is expected.

We also compare two regions in the spectrum in which the spectroscopic signatures of NH₃ are visually prominent. These are features at about 3 µm with NIRCam F322W2 and about 10 µm with MIRI. These two features are parts of regions in which NH₃ contributes >50% of
the total cross-section (see shaded purple regions in Extended Data Fig. 8). The region between 2.9 and 3.1 µm covered by NIRCam F322W2 has an order of magnitude greater density in the increased predictive performance of NH₃ than the 8.6–12 µm region covered by MIRI. Our analysis finds that the detection of NH₃ is more strongly driven by the NIRCam observations than the MIRI observations.

NIRCam F322W2 provides the information necessary to confirm the tentative detection (approximately 2–3σ) of NH₃ suggested by ref. 16 using MIRI data alone. Our LOO-CV analysis suggests that definitive detections of NH₃ require resolving the strong spectroscopic feature at about 3.0 µm and highlight an important advantage of NIRCam over other instruments aboard the JWST that do not have the necessary wavelength coverage or throughput. Although NH₃ has a strong spectral feature at roughly 10 µm, and the cross-section of the gas is dominant at >8.6 µm, the presence of strong cloud resonance features in the infrared obfuscates the observational ability to strongly detect the gas.

**Stellar SED fitting and absolute radii.** We used a set of catalogue magnitudes for the WASP-107 system to fit a model spectral energy distribution (SED) to the stellar emission. In conjunction with the Gaia¹¹² parallax for the system, this allowed us to estimate a stellar radius for the star WASP-107.

For our SED fits, we used catalogue 2MASS JHK¹¹³ and AllWISE W1 and W2 (ref. 114) magnitudes. Our SED model used four different physical parameters: the stellar effective temperature, the stellar radius, the amount of visual extinction to WASP-107 and the parallax of the system. We assumed $\log g = 4.5$ and $[\text{Fe/H}] = 0.0$ for the star WASP-107, which are both within 0.1 dex of the surface gravity and metallicity measured previously⁵. In practice, the precise values of logg and [Fe/H] do not greatly affect the results from the SED fit¹¹⁵.

We imposed a Gaussian prior on $T_{\text{eff}}$ based on previous spectroscopic measurements⁵ of $T_{\text{eff}} = 4,425 \pm 70$ K, with the associated 10σ uncertainties as the prior width. Similarly, we imposed a Gaussian prior on the parallax to the system using the Gaia DR3 (ref. 112) parallax of $\pi = 15.528 \pm 0.026$ mas. For the amount of visual extinction to WASP-107, we imposed a prior of $A_V = 0.03 \pm 0.01$, which comes from measurements¹¹⁶ of the excess reddening towards WASP-107 of $E(B-V) = 0.028 \pm 0.01$, and assuming $R_V = 3.1$.

To model the SED, we used BHAC15 spectra¹¹⁷ for the stellar SED. We computed a grid of surface luminosity magnitudes, corresponding to the bandpasses of the catalogue magnitudes, for a range of $T_{\text{eff}}$ values. Because the BHAC15 models step by 200 K in $T_{\text{eff}}$, we used cubic spline interpolation to estimate model magnitudes between the points provided by the model atmospheres. We then scaled the interpolated surface magnitudes for the star by $R_i/d$–in which *d* is the distance to the star—to determine the apparent bolometric flux of the SED at Earth. We then applied a simple $R = 3.1$ extinction law scaled from the value of $A_V$ to determine the extincted bolometric flux of the SED model.

This stellar SED fitting allowed us to measure the radius of the star WASP-107 to be $R_⊙ = 0.67 \pm 0.01 R_V$. We then used the mean value of $R_p/R_⊙$ as measured in our joint broadband transit fits to estimate the planetary radius of WASP-107b to be $R_p = 0.939 \pm 0.019 R_⊙$.

## Data availability

The NIRCam data used in this paper are associated with JWST GTO programme 1185 (P.I. Greene; observations 8 and 9) and will be publicly available from the Mikulski Archive for Space Telescopes (MAST; https://mast.stsci.edu) at the end of their 1-year exclusive access period. The MIRI data used in this paper are from JWST GTO programme 1280 (P.I. Lagage; observation 1) and will also be publicly available on the MAST at the end of their proprietary period. Further intermediate results from this work are archived on Zenodo at https://doi.org/10.5281/zenodo.10780448 (ref. 46).

## Code availability

We used the following codes to reduce and fit the JWST data: STScI's JWST calibration pipeline⁴⁴, Eureka!¹⁴, tshirt¹⁵, starry⁶³, PyMC3 (ref. 66), numpy¹¹⁸, astropy¹¹⁹,¹²⁰, scipy¹²¹ and matplotlib¹²².

44. Bushouse, H. et al. JWST calibration pipeline. Zenodo https://doi.org/10.5281/zenodo.7829392 (2023).

45. Ahrer, E.-M. et al. Early Release Science of the exoplanet WASP-39b with JWST NIRCam. Nature **614**, 653–658 (2023).

46. Welbanks, L. et al. Products and data for “A high internal heat flux and large core for the inflated warm Neptune WASP-107b”. Zenodo https://doi.org/10.5281/zenodo.10780448 (2024).

47. Horne, K. An optimal extraction algorithm for CCD spectroscopy. Publ. Astron. Soc. Pac. **98**, 609–617 (1986).

48. Schlawin, E. et al. JWST noise floor. I. Random error sources in JWST NIRCam time series. Astron. J. **160**, 231 (2020).

49. Schlawin, E. et al. JWST NIRCam defocused imaging: photometric stability performance and how it can sense mirror tilts. Publ. Astron. Soc. Pac. **135**, 018001 (2023).

50. Beatty, T. G. et al. Evidence for atmospheric cold-trop processes in the noninverted emission spectrum of Kepler-13Ab using HST/WFC3. Astron. J. **154**, 158 (2017).

51. Bell, T. J. et al. Nightside clouds and disequilibrium chemistry on the hot Jupiter WASP-43b. Nat. Astron. https://doi.org/10.1038/s41550-024-02230-x (2024).

52. Bouwman, J. et al. Spectroscopic time series performance of the Mid-infrared Instrument on the JWST. Publ. Astron. Soc. Pac. **135**, 038002 (2023).

53. Lomb, N. R. Least-squares frequency analysis of unequally spaced data. Astrophys. Space Sci. **39**, 447–462 (1976).

54. Scargle, J. D. Studies in astronomical time series analysis. II. Statistical aspects of spectral analysis of unevenly spaced data. Astrophys. J. **263**, 835–853 (1982).

55. Kendrew, S. et al. The Mid-Infrared Instrument for the James Webb Space Telescope, IV: the low-resolution spectrometer. Publ. Astron. Soc. Pac. **127**, 623 (2015).

56. Bell, T. J. et al. A first look at the JWST MIRI/LRS phase curve of WASP-43b. Zenodo https://zenodo.org/doi/10.5281/zenodo.7547783 (2023).

57. Clemens, J. C., Crain, J. A. & Anderson, R. The Goodman spectrograph. Proc. SPIE **5492**, 331–340 (2004).

58. Vogt, S. S. et al. HIRES: the high-resolution echelle spectrometer on the Keck 10-m Telescope. Proc. SPIE **2198**, 362–375 (1994).

59. Collins, K. A., Kielkopf, J. F., Stassun, K. G. & Hessman, F. V. Astrolmage: image processing and photometric extraction for ultra-precise astronomical light curves. Astron. J. **153**, 77 (2017).

60. Beatty, T. G. et al. A significant overluminosity in the transiting brown dwarf CWW 89Ab. Astron. J. **156**, 168 (2018).

61. Kreidberg, L. batman: BAsic Transit Model cAlculationN in Python. Publ. Astron. Soc. Pac. **127**, 1161 (2015).

62. Foreman-Mackey, D., Hogg, D. W., Lang, D. & Goodman, J. emcee: the MCMC hammer. Publ. Astron. Soc. Pac. **125**, 306 (2013).

63. Luger, R. et al. starry: analytic occultation light curves. Astron. J. **157**, 64 (2019).

64. Kurucz, R. L. ATLAS9 stellar atmosphere programs and 2km/s grid. Kurucz CD-Rom 13 (Smithsonian Astrophysical Observatory, 1993).

65. Bourque, M. et al. The Exoplanet Characterization Toolkit (ExoCTK). Zenodo https://doi.org/10.5281/zenodo.4556063 (2021).

66. Salvatier, J., Wiecki, T. V. & Fonnesbeck, C. Probabilistic programming in Python using PyMC3. PeerJ Comput. Sci. **2**, e55 (2016).

67. Gelman, A. & Rubin, D. B. Inference from iterative simulation using multiple sequences. Stat. Sci. **7**, 457–472 (1992).

68. Allan, D. W. Statistics of atomic frequency standards. Proc. IEEE **54**, 221–230 (1966).

69. Winn, J. N. et al. The transit light curve project. IX. Evidence for a smaller radius of the exoplanet XO-3b. Astrophys. J. **683**, 1076–1084 (2008).

70. Kipping, D. M. Efficient, uninformative sampling of limb darkening coefficients for two-parameter laws. Mon. Not. R. Astron. Soc. **435**, 2152–2160 (2013).

71. Magic, Z., Chiavassa, A., Collet, R. & Asplund, M. The Stagger-grid: a grid of 3D stellar atmosphere models. IV. Limb darkening coefficients. Astron. Astrophys. **573**, A90 (2015).

72. Grant, D. & Wakeford, H. R. Exo-TiC/ExoTiC-LD: ExoTiC-LD v3.0.0. Zenodo https://doi.org/10.5281/zenodo.7437681 (2022).

73. Piskorz, D. et al. Ground- and space-based detection of the thermal emission spectrum of the transiting hot Jupiter KELT-2Ab. Astron. J. **156**, 133 (2018).

74. Mansfield, M. et al. A unique hot Jupiter spectral sequence with evidence for compositional diversity. Nat. Astron. **5**, 1224–1232 (2021).

75. Tsai, S.-M. et al. VULCAN: an open-source, validated chemical kinetics Python code for exoplanetary atmospheres. Astrophys. J. Suppl. Ser. **228**, 20 (2017).

76. Husser, T. O. et al. A new extensive library of PHOENIX stellar atmospheres and synthetic spectra. Astron. Astrophys. **553**, A6 (2013).

77. France, K. et al. The MUSCLES Treasury Survey. I. Motivation and overview. Astrophys. J. **820**, 89 (2016).

78. Youngblood, A. et al. The MUSCLES Treasury Survey. II. Intrinsic LYα and extreme ultraviolet spectra of K and M dwarfs with exoplanets. Astrophys. J. **824**, 101 (2016).

79. Loyd, R. O. P. et al. The MUSCLES Treasury Survey. III. X-ray to infrared spectra of 11 M and K stars hosting planets. Astrophys. J. **824**, 102 (2016).

80. Welbanks, L. & Madhusudhan, N. On degeneracies in retrievals of exoplanetary transmission spectra. Astron. J. **157**, 206 (2019).

81. Underwood, D. S. et al. ExoMol molecular line lists – XIV. The rotation–vibration spectrum of hot SO₂. Mon. Not. R. Astron. Soc. **459**, 3890–3899 (2016).

82. Benneke, B. & Seager, S. How to distinguish between cloudy mini-Neptunes and water/volatile-dominated super-Earths. Astrophys. J. **778**, 153 (2013).
83. Gao, P., Wakeford, H. R., Moran, S. E. & Parmentier, V. Aerosols in exoplanet atmospheres. *J. Geophys. Res. Planets* **126**, e06655 (2021).

84. Line, M. R. & Parmentier, V. The influence of nonuniform cloud cover on transit transmission spectra. *Astrophys. J.* **820**, 78 (2016).

85. Welbanks, L. & Madhusudhan, N. On atmospheric retrievals of exoplanets with inhomogeneous terminators. *Astrophys. J.* **933**, 79 (2022).

86. Barstow, J. K. Unveiling cloudy exoplanets: the influence of cloud model choices on retrieval solutions. *Mon. Not. R. Astron. Soc.* **497**, 4183–4195 (2020).

87. Tsiaras, A. et al. A population study of gaseous exoplanets. *Astron. J.* **155**, 156 (2018).

88. Fisher, C. & Heng, K. Retrieval analysis of 38 WFC3 transmission spectra and resolution of the normalization degeneracy. *Mon. Not. R. Astron. Soc.* **481**, 4698–4727 (2018).

89. Benneke, B. et al. A sub-Neptune exoplanet with a low-metallicity methane-depleted atmosphere and Mie-scattering clouds. *Nat. Astron.* **3**, 813–821 (2019).

90. Pinhas, A., Rackham, B. V., Madhusudhan, N. & Apai, D. Retrieval of planetary and stellar properties in transmission spectroscopy with AURA. *Mon. Not. R. Astron. Soc.* **480**, 5314–5331 (2018).

91. Welbanks, L. & Madhusudhan, N. Aurora: a generalized retrieval framework for exoplanetary transmission spectra. *Astrophys. J.* **913**, 114 (2021).

92. Grant, D. et al. JWST-TST DREAMS: quartz clouds in the atmosphere of WASP-17b. *Astrophys. J. Lett.* **956**, L29 (2023).

93. Ackerman, A. S. & Marley, M. S. Precipitating condensation clouds in substellar atmospheres. *Astrophys. J.* **556**, 872–884 (2001).

94. Mollière, P. et al. petitRADTRANS: A Python radiative transfer package for exoplanet characterization and retrieval. *Astron. Astrophys.* **627**, A67 (2019).

95. Wakeford, H. R. & Sing, D. K. Transmission spectral properties of clouds for hot Jupiter exoplanets. *Astron. Astrophys.* **573**, A122 (2015).

96. Mai, C. & Line, M. R. Exploring exoplanet cloud assumptions in JWST transmission spectra. *Astrophys. J.* **883**, 144 (2019).

97. Line, M. R., Knutson, H., Deming, D., Wilkins, A. & Desert, J.-M. A near-infrared transmission spectrum for the warm Saturn HAT-P-12b. *Astrophys. J.* **778**, 183 (2013).

98. Madhusudhan, N. & Seager, S. A temperature and abundance retrieval method for exoplanet atmospheres. *Astrophys. J.* **707**, 24–39 (2009).

99. Skilling, J. in *AIP Conference Proceedings*, Vol. 735 (eds Fischer, R., Preuss, R. & Toussaint, U. V.) 395–405 (American Institute of Physics, 2004).

100. Feroz, F., Hobson, M. P. & Bridges, M. MULTINEST: an efficient and robust Bayesian inference tool for cosmology and particle physics. *Mon. Not. R. Astron. Soc.* **398**, 1601–1614 (2009).

101. Buchner, J. et al. X-ray spectral modelling of the AGN obscuring region in the CDFS: Bayesian model selection and catalogue. *Astron. Astrophys.* **564**, A125 (2014).

102. Feroz, F., Hobson, M. P., Cameron, E. & Pettitt, A. N. Importance nested sampling and the MultiNest algorithm. *Open J. Astrophys.* **2**, 10 (2019).

103. Rothman, L. S. et al. HITEMP, the high-temperature molecular spectroscopic database. *J. Quant. Spectrosc. Radiat. Transf.* **111**, 2139–2150 (2010).

104. Yurchenko, S. N. & Tennyson, J. ExoMol line lists – IV. The rotation–vibration spectrum of methane up to 1500 K. *Mon. Not. R. Astron. Soc.* **440**, 1649–1661 (2014).

105. Yurchenko, S. N., Barber, R. J. & Tennyson, J. A variationally computed line list for hot NH₃. *Mon. Not. R. Astron. Soc.* **413**, 1828–1834 (2011).

106. Welbanks, L., McGill, P., Line, M. & Madhusudhan, N. On the application of Bayesian leave-one-out cross-validation to exoplanet atmospheric analysis. *Astron. J.* **165**, 112 (2023).

107. Vehtari, A. & Ojanen, J. A survey of Bayesian predictive methods for model assessment, selection and comparison. *Stat. Surv.* **6**, 142–228 (2012).

108. Vehtari, A., Gelman, A. & Gabry, J. Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. *Stat. Comput.* **27**, 1413–1432 (2017).

109. Vehtari, A., Simpson, D., Gelman, A., Yao, Y. & Gabry, J. Pareto smoothed importance sampling. *J. Mach. Learn. Res.* **25**, 1–58 (2024).

110. McGill, P. et al. First semi-empirical test of the white dwarf mass–radius relationship using a single white dwarf via astrometric microlensing. *Mon. Not. R. Astron. Soc.* **520**, 259–280 (2023).

111. Challener, R. C., Welbanks, L. & McGill, P. Bringing 2D eclipse mapping out of the shadows with leave-one-out cross validation. *Astron. J.* **166**, 251 (2023).

112. Gaia Collaboration. Gaia Data Release 3. Summary of the content and survey properties. *Astron. Astrophys.* **674**, A1 (2023).

113. Skrutskie, M. F. et al. The Two Micron All Sky Survey (2MASS). *Astron. J.* **131**, 1163–1183 (2006).

114. Cutri, R. M. et al. VizieR Online Data Catalog: AllWISE Data Release (Cutri+ 2013). VizieR On-line Data Catalog: II/328 (IPAC/Caltech, 2014).

115. Stevens, D. J., Gaudí, B. S. & Stassun, K. G. Measuring model-independent masses and radii of single-lined eclipsing binaries: analytic precision estimates. *Astrophys. J.* **862**, 53 (2018).

116. Schlafly, E. F. & Finkbeiner, D. P. Measuring reddening with Sloan Digital Sky Survey stellar spectra and recalibrating SFD. *Astrophys. J.* **737**, 103 (2011).

117. Baraffe, I., Homeier, D., Allard, F. & Chabrier, G. New evolutionary models for pre-main sequence and main sequence low-mass stars down to the hydrogen-burning limit. *Astron. Astrophys.* **577**, A42 (2015).

118. Harris, C. R. et al. Array programming with NumPy. *Nature* **585**, 357–362 (2020).

119. Astropy Collaboration. Astropy: a community Python package for astronomy. *Astron. Astrophys.* **558**, A33 (2013).

120. Astropy Collaboration. The Astropy Project: building an open-science project and status of the v2.0 core package. *Astron. J.* **156**, 123 (2018).

121. Vintanen, P. et al. SciPy 1.0: fundamental algorithms for scientific computing in Python. *Nat. Methods* **17**, 261–272 (2020).

122. Hunter, J. D. Matplotlib: a 2D graphics environment. *Comput. Sci. Eng.* **9**, 90–95 (2007).

**Acknowledgements** L.W. acknowledges support for this work provided by NASA through the NASA Hubble Fellowship grant #HST-HF2-51496.001-A awarded by the Space Telescope Science Institute (STScI), which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. T.P.G. and T.J.B. acknowledge support from NASA JWST WBSs 411672.07.04.01.02 and 411672.07.05.05.03.02. P.M. acknowledges that this work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under contract DE-AC52-07NA27344. The document number is LLNL-JRNL-859050. M.R.L. acknowledges NASA XRP award 80NSSC19K0446 and STScI grant HST-AR-16139. M.R.L. and L.W. acknowledge Research Computing at Arizona State University for providing high-performance computing and storage resources that have substantially contributed to the research results reported in this manuscript. L.W. thanks M. Min for meaningful conversations. K.O. acknowledges support from JSPS KAKENHI grant number JP23K19072. M.M. acknowledges funding from NASA Goddard Spaceflight Center through NASA contract NAS5-02105. A.D. and P.-O.L. acknowledge funding support from CNES. We thank M. Rieke for allocating the NIRCam time for this programme.

**Author contributions** L.W. led the modelling analysis effort, performed the atmospheric modelling using free retrievals and 1D-RCPE models, contributed to the theoretical analysis/ interpretation of the observations, led the cross-validation analysis and led the writing of the manuscript. T.J.B. contributed the fiducial Eureka! analyses of the NIRCam and MIRI observations, verified the observing parameters of the NIRCam observations and contributed to the text. T.G.B. contributed the Pegasus analyses of the HST and NIRCam observations, verified the observing parameters of the NIRCam observations, performed the SED fitting and radius estimation, performed the tidal heating analysis and contributed to the text. M.R.L. performed the 1D-RCPE simulations, wrote the introductory text and methods and helped sculpt the direction of the manuscript. K.O. provided comments on the manuscript and helped with interpreting the NH₃ detection. J.J.F. contributed to the planet structure models, contributed to the interpretation of the internal temperature of the planet and contributed to the text. E.S. simulated the planet and retrievals before launch, help planned the observation specifications, carried out the tshirt reduction and light-curve fitting and contributed to the text. T.P.G. selected the planet for observation, designed the observational programme, directed some of the analysis and commented on the manuscript. E.R. provided feedback on the modelling interpretation and provided comments on the manuscript. P.M. contributed the cross-validation analysis and to the text. V.P. provided feedback throughout the project and helped with the vertical mixing discussion. Y.T. contributed to the planet structure models. M.M. contributed to deriving the orbital parameters of the planet and feedback on the text. I.E., S.M., L.S.W. and K.E.A. provided useful discussions throughout the manuscript process. P.-O.L. designed the MIRI observational programme and provided its data. A.D. contributed to the MIRI data analysis.

**Competing interests** The authors declare no competing interests.

**Additional information**

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41586-024-07514-w.

**Correspondence and requests for materials** should be addressed to Luis Welbanks.

**Peer review information** Nature thanks the anonymous reviewers for their contribution to the peer review of this work. Peer reviewer reports are available.

**Reprints and permissions information** is available at http://www.nature.com/reprints.
<figure><img src="image_7.png" /></figure>

**Extended Data Fig. 1** | The 390-Hz noise signal of MIRI/LRS visualized.
a, The background pixel values for part of the second integration of the MIRI/LRS observations of WASP-107b are shown by black points, whereas our fitted 390-Hz noise signal is shown with a red line. b, The Lomb–Scargle periodogram of the pixels in the second group of the second integration before the 390-Hz noise-removal and GLBS steps are applied are shown in red and the periodogram

<figure><img src="image_8.png" /></figure>

after these two steps are applied is shown in black. The three harmonics used in our cleaning procedure (1st, 2nd and 4th) are indicated with pale blue vertical lines and the skipped 3rd harmonic is shown with a pale orange vertical line. The clear impact of the 390-Hz noise-removal and GLBS steps is seen in the eliminated spikes in the Lomb–Scargle periodogram.
<figure><img src="image_9.png" /></figure>

**Extended Data Fig. 2 | Comparison of different NIRCam reductions and limb-darkening choices.** a, The NIRCam F322W2 and F444W transmission spectra as reduced by the Eureka!, Pegasus and tshirt pipelines are shown in different colours after being smoothed by a four-point-wide boxcar filter. Both the Eureka! and Pegasus reductions fixed limb-darkening coefficients to ATLAS model predictions, whereas tshirt freely fit quadratic coefficients. At wavelengths longer than about 2.9 µm, the F322W2 spectrum of tshirt

<figure><img src="image_10.png" /></figure>

consistently falls below the spectra of Eureka! and Pegasus; this difference is caused by a combination of tshirt's ROEBA method (which results in a roughly constant offset across F322W2) and tshirt's free limb darkening (which results in a small slope across F322W2). b, A similar plot, demonstrating the minimal impacts of varying limb-darkening choices on our final spectra produced by Eureka!
L168-9b MIRI/LRS commissioning observations

<figure><img src="image_11.png" /></figure>

**Extended Data Fig. 3 | A demonstration of the underestimation of error bars in the MIRI/LRS observations of L168-9b and the impact of the 390-Hz noise-removal step.** The median fitted transit depth uncertainties decrease with increasing spectral bin width, as would be expected for white noise (blue symbols). Meanwhile, the standard deviation of the L168-9b transmission

spectrum is around two times the expected noise level in the spectrally unbinned data. The ratio of the standard deviation of the spectrum to the median fitted uncertainty decreases with increasing bin size, with the exception of a peak around 0.25 µm.
<figure><img src="image_12.png" /><figcaption>**Extended Data Fig. 4** | A comparison of different reductions of the WASP-107b MIRI/LRS spectra. a, A demonstration of the impact of the 390-Hz noise-removal and GLBS steps on our reduction of the WASP-107b MIRI/LRS data. For wavelengths >10 µm, the GLBS step slightly reduces the final transmission, while also turning off the 390-Hz noise-removal step seems to have no further impact at any wavelength. Error bars show 1σ uncertainties. b, A comparison of our fiducial reduction with the three reductions presented in ref. 16, with 1σ uncertainties. Although there are some differences between the four different</figcaption></figure>

reductions, the overall shape of the transmission spectrum is robust to
differences in reduction choices. c, The per-point differences between different
pairs of reductions, which are summarized with a histogram in panel d. With
the exception of a single TEATRO point, our reduction agrees with the European
Consortium's three reductions at ≲3σ. However, there is some structure to the
differences between our reduction and those of the European Consortium,
with our reduction giving larger transit depths on average than the European
Consortium's from around 8–10 µm.
<figure><img src="image_13.png" /></figure>

Extended Data Fig. 5 | Illustration of the effects of the key physical processes on the atmospheric structure and resultant spectra. Panels a and b show the impact of disequilibrium chemical processes (vertical mixing and photochemistry) relative to thermochemical equilibrium. Panels c and d show the influence of changing internal temperature. Panels e and f show the impact

of changing eddy diffusion strength. We only show the vertical abundance profiles (panels a, c and e) for observable gases that are influenced by these effects (for instance, water is not shown as it not affected). The observed transit spectrum is shown with grey points with 1σ error bars in panels b, d and f.
<figure><img src="image_14.png" /></figure>

**Extended Data Fig. 6 | Effect of self-consistent microphysical eddysed clouds on the spectrum.** Panel a shows the condensate vertical distributions for the major condensate species. The cloud bases, in which the condensates

<figure><img src="image_15.png" /></figure>

first forms, all occur at or below the 100-mbar level, with the silicate clouds forming at or below the 1-bar level. Panel b shows the resulting spectrum (with 1σ error bars), as well as the contribution from each cloud species.
<figure><img src="image_16.png" /></figure>

**Extended Data Fig. 7** | Comparison of retrieved pressure-temperature structure. The retrieved vertical temperature structures (median, 1σ and 2σ) from the free retrieval (Aurora, blue) and the 1D-RCPE model (red) are generally in good agreement within their 68% confidence intervals. The dotted line shows

the equilibrium temperature of the planet. The observations generally examine pressures lower than about 100 mbar and as low as about several 10⁻⁵ bar. We include in black the T_int = 200 K profile from Extended Data Fig. 5 as an example of lower internal temperatures.
<figure><img src="image_17.png" /></figure>

**Extended Data Fig. 8** | The detection of NH₃ in the transmission spectrum of WASP-107b is driven by NIRCAM F322W2 observations. The data are colour-coded by the point-wise difference in expected log posterior predictive density between the reference 1D-RCPE model and the model without NH₃ absorption. The best-fit 1D-RCPE model (*R* = 300) is shown in grey. Redder data points (larger positive Δelpd score) are better explained by the reference model including NH₃ absorption. Although points with positive Δelpd are present throughout the entire transmission spectrum, indicating that NH₃ improves

the model performance at all wavelengths, the highest scoring points are localized at roughly 3 µm, at which strong NH₃ absorption is visible. The purple-shaded regions show areas in which the cross-section of NH₃ contributes to 50% or more of the cross-sections detected in WASP-107b, and with visually prominent features. Three data points in the MIRI observations (7.017, 7.318 and 7.468 µm) are not included in the analysis because their associated Pareto *k* value exceeded 0.7. Data points are shown with 1σ error bars.
## Extended Data Table 1 | 1D-RCPE retrieved atmospheric properties

<table><thead><tr><th></th><th>Parameter</th><th>Value</th><th>Prior</th></tr></thead><tbody><tr><td rowspan="7">Atmospheric State</td><td>[M/H]</td><td>1.09<sup>+0.17</sup><sub>-0.07</sub></td><td>(0.5, 2.0)</td></tr><tr><td>C/O</td><td>0.33<sup>+0.06</sup><sub>-0.05</sub></td><td>(0.05, 0.60)</td></tr><tr><td>T<sub>int</sub></td><td>493<sup>+41</sup><sub>-56</sub></td><td>(200, 550)</td></tr><tr><td>log<sub>10</sub>(K<sub>zz</sub>)</td><td>8.62<sup>+0.42</sup><sub>-0.22</sub></td><td>(7, 10)</td></tr><tr><td>log<sub>10</sub>(&alpha;NH<sub>3</sub>)</td><td>1.46<sup>+0.17</sup><sub>-0.19</sub></td><td>(-4, 4)</td></tr><tr><td>&times; R<sub>p</sub></td><td>1.03<sup>+0.03</sup><sub>-0.03</sub></td><td>(0.5, 1.5)</td></tr><tr><td>log<sub>10</sub>(P<sub>ref.</sub>) (bar)</td><td>-3.76<sup>+1.5</sup><sub>-1.4</sub></td><td>(-6.0, 1.5)</td></tr><tr><td rowspan="8">Cloud/Haze</td><td>log<sub>10</sub>(K<sub>gray clouds</sub>)</td><td>-38.0<sup>+4.8</sup><sub>-4.6</sub></td><td>(-45, -20)</td></tr><tr><td>log<sub>10</sub>(a)</td><td>2.3<sup>+0.2</sup><sub>-0.2</sub></td><td>(-14, 14)</td></tr><tr><td>&gamma;</td><td>-1.3<sup>+0.1</sup><sub>-0.1</sub></td><td>(-20, 2)</td></tr><tr><td>log<sub>10</sub>(K<sub>cond</sub>)</td><td>-28.0<sup>+0.2</sup><sub>-0.2</sub></td><td>(-45, -20)</td></tr><tr><td>&lambda;<sub>0</sub></td><td>8.2<sup>+0.7</sup><sub>-0.6</sub></td><td>(5, 13)</td></tr><tr><td>log<sub>10</sub>(&omega;)</td><td>0.07<sup>+0.08</sup><sub>-0.05</sub></td><td>(-1, 1)</td></tr><tr><td>&xi;</td><td>0.4<sup>+1.0</sup><sub>-0.9</sub></td><td>(-10, 10)</td></tr><tr><td>&phi;<sub>clouds and hazes</sub></td><td>0.9<sup>+0.1</sup><sub>-0.1</sub></td><td>(0, 1)</td></tr><tr><td></td><td>&Delta;<sub>MIRI</sub> (ppm)</td><td>282<sup>+49</sup><sub>-41</sub></td><td>(-50, 500)</td></tr></tbody></table>

Retrieved parameters (median and 1σ uncertainties) and their prior ranges are included. log₁₀(αNH₃) is the enhancement to the NH₃ abundance profile and ΔMIRI is the offset between the MIRI observations and NIRCAM F444W.
Extended Data Table 2 | Free retrieval retrieved atmospheric properties

<table><thead><tr><th></th><th>Parameter</th><th>CHIMERA</th><th>Aurora</th><th>Prior</th></tr></thead><tbody><tr><td rowspan="10">Chemical Species</td><td>log<sub>10</sub>(X<sub>H2O</sub>)</td><td>-2.1<sup>+0.2</sup><sub>-0.3</sub></td><td>-2.6<sup>+0.7</sup><sub>-0.6</sub></td><td>(-12.0, -0.3)</td></tr><tr><td>log<sub>10</sub>(X<sub>CH4</sub>)</td><td>-5.8<sup>+0.2</sup><sub>-0.2</sub></td><td>-6.1<sup>+0.5</sup><sub>-0.4</sub></td><td>(-12.0, -0.3)</td></tr><tr><td>log<sub>10</sub>(X<sub>NH3</sub>)</td><td>-5.0<sup>+0.2</sup><sub>-0.2</sub></td><td>-5.1<sup>+0.5</sup><sub>-0.5</sub></td><td>(-12.0, -0.3)</td></tr><tr><td>log<sub>10</sub>(X<sub>CO</sub>)</td><td>-1.9<sup>+0.2</sup><sub>-0.2</sub></td><td>-3.0<sup>+0.7</sup><sub>-0.7</sub></td><td>(-12.0, -0.3)</td></tr><tr><td>log<sub>10</sub>(X<sub>CO2</sub>)</td><td>-3.9<sup>+0.3</sup><sub>-0.3</sub></td><td>-4.4<sup>+0.8</sup><sub>-0.7</sub></td><td>(-12.0, -0.3)</td></tr><tr><td>log<sub>10</sub>(X<sub>SO2</sub>)</td><td>-5.2<sup>+0.2</sup><sub>-0.2</sub></td><td>-5.7<sup>+0.4</sup><sub>-0.4</sub></td><td>(-12.0, -0.3)</td></tr><tr><td>log<sub>10</sub>(X<sub>H2S</sub>)</td><td>-8.5<sup>+2.2</sup><sub>-2.2</sub></td><td>-8.6<sup>+2.1</sup><sub>-2.2</sub></td><td>(-12.0, -0.3)</td></tr><tr><td>T<sub>0</sub> (K)</td><td>525<sup>+26</sup><sub>-31</sub></td><td>447<sup>+42</sup><sub>-30</sub></td><td>(400, 900)</td></tr><tr><td>α<sub>1</sub> (K<sup>-1/2</sup>)</td><td>1.6<sup>+0.3</sup><sub>-0.3</sub></td><td>1.6<sup>+0.3</sup><sub>-0.3</sub></td><td>(0.02, 2.0)</td></tr><tr><td>α<sub>2</sub> (K<sup>-1/2</sup>)</td><td>1.2<sup>+0.6</sup><sub>-0.6</sub></td><td>1.2<sup>+0.5</sup><sub>-0.6</sub></td><td>(0.02, 2.0)</td></tr><tr><td rowspan="4">P-T</td><td>log<sub>10</sub>(P<sub>1</sub>) (bar)</td><td>-2.4<sup>+1.6</sup><sub>-1.9</sub></td><td>-2.0<sup>+1.7</sup><sub>-1.9</sub></td><td>(-9.0, 2.0)</td></tr><tr><td>log<sub>10</sub>(P<sub>2</sub>) (bar)</td><td>-5.8<sup>+2.0</sup><sub>-1.8</sub></td><td>-5.8<sup>+2.5</sup><sub>-2.1</sub></td><td>(-9.0, 2.0)</td></tr><tr><td>log<sub>10</sub>(P<sub>3</sub>) (bar)</td><td>-0.1<sup>+0.9</sup><sub>-1.0</sub></td><td>0.4<sup>+1.1</sup><sub>-1.3</sub></td><td>(-2.0, 2.0)</td></tr><tr><td>log<sub>10</sub>(P<sub>rel</sub>) (bar)</td><td>-4.9<sup>+2.4</sup><sub>-2.0</sub></td><td>-7.3<sup>+1.3</sup><sub>-1.0</sub></td><td>(-9.0, 2.0)</td></tr><tr><td></td><td>R<sub>p</sub> (R<sub>Jup</sub>)</td><td>0.98<sup>+0.04</sup><sub>-0.05</sub></td><td>1.01<sup>+0.02</sup><sub>-0.03</sub></td><td>(0.66, 1.22)</td></tr><tr><td rowspan="10">Cloud/Haze</td><td>log<sub>10</sub>(a)</td><td>2.5<sup>+0.2</sup><sub>-0.2</sub></td><td>1.7<sup>+0.5</sup><sub>-0.4</sub></td><td>(-4, 10)</td></tr><tr><td>γ</td><td>-1.0<sup>+0.1</sup><sub>-0.1</sub></td><td>-1.2<sup>+0.1</sup><sub>-0.2</sub></td><td>(-20, 2)</td></tr><tr><td>log<sub>10</sub>(K<sub>cond</sub>)</td><td>-27.9<sup>+0.2</sup><sub>-0.2</sub></td><td>-28.4<sup>+0.5</sup><sub>-0.3</sub></td><td>(-50, -20)</td></tr><tr><td>λ<sub>0</sub></td><td>7.7<sup>+0.5</sup><sub>-0.4</sub></td><td>7.4<sup>+0.5</sup><sub>-0.3</sub></td><td>(5.5, 12.0)</td></tr><tr><td>log<sub>10</sub>(ω)</td><td>0.1<sup>+0.1</sup><sub>-0.1</sub></td><td>0.2<sup>+0.1</sup><sub>-0.1</sub></td><td>(-1, 1)</td></tr><tr><td>ξ</td><td>1.1<sup>+0.9</sup><sub>-0.7</sub></td><td>2.3<sup>+1.5</sup><sub>-1.2</sub></td><td>(-10, 10)</td></tr><tr><td>log<sub>10</sub>(P<sub>cloud</sub>) (bar)</td><td>N/A</td><td>-0.1<sup>+1.3</sup><sub>-1.3</sub></td><td>(-9.0, 2.0)</td></tr><tr><td>log<sub>10</sub>(κ<sub>cloud</sub>)</td><td>-42.8<sup>+8.1</sup><sub>-7.4</sub></td><td>N/A</td><td>(-55, -20)</td></tr><tr><td>φ<sub>clouds and hazes</sub></td><td>0.9<sup>+0.1</sup><sub>-0.1</sub></td><td>0.9<sup>+0.1</sup><sub>-0.1</sub></td><td>(0, 1)</td></tr><tr><td>Δ<sub>MIRI</sub> (ppm)</td><td>253<sup>+41</sup><sub>-39</sub></td><td>285<sup>+37</sup><sub>-36</sub></td><td>(-50, 500)</td></tr></tbody></table>

Retrieved parameters (median and 1σ uncertainties) and their prior ranges are included.