A High Internal Heat Flux and Large Core in a Warm Neptune Exoplanet (2024)

Eureka!

Our Eureka! reduction used version 0.9 of the Eureka! pipeline[14], CRDS version 11.17.0 and context ‘1093’ (for F322W2) or ‘1097’ (for F444W), and jwst package version 1.10.2 (ref.\citeAppjwst_v1.10.2). Our reduction methods closely follow those used in previous Eureka! NIRCam spectroscopy analyses[11]^,\citeAppahrer2022nircam. The Eureka! Control Files and Eureka! Parameter Files we used are available for download from Zenodo (https://doi.org/10.5281/zenodo.10780448; ref.\citeAppwelbanks2024_wasp107b_zenodo), 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$\sigma$ from the default of 4$\sigma$ (to avoid excessive false positives) and turning off the photom step in Stage 2 (as flux-calibrated data is not desired for time-series observations). In Eureka!’s Stage 3, we first did two iterations of 5$\sigma$ clipping along the time axis on the background pixels, and we also masked pixels marked as ‘DO_NOT_USE’ in the data quality (DQ) 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 centered 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 center of the star’s trace and masking any 7$\sigma$ outliers in the column. The spatial position and 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 lightcurves in Eureka!’s Stage 5. We then performed optimal spectral extraction\citeApphorne1986optspec using only the pixels within 5 pixels of the center 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$\sigma$ outliers along the time axis and then smoothed along the spectral direction using a boxcar filter with a width of 13 pixels. Pixels which differed by more than 10$\sigma$ compared to 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 (likely due to 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 $\mu$m wide bins for both NIRCam filters (100 bins spanning 2.45–3.95 $\mu$m for F322W2 and 72 bins spanning 3.89–4.97 $\mu$m for F444W). Finally, to remove any remaining cosmic ray effects, we clipped any 4$\sigma$ outliers in the spectrally binned lightcurves compared to 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 prep.), and 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$\times$2,048 rateints images. We masked out image rows 5 to 75 to not self-subtract light from the WASP-107 system, and then we fit 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 re-ran the reference pixel correction using hxrg-ref-pixel (https://github.com/JarronL/hxrg_ref_pixels). This appeared to correct the issue.

We then performed a spectrophotometric extraction on the background-subtracted images using optimal extraction techniques\citeApphorne1986optspec. 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\citeApphorne1986optspec. We then fit a 4th-order polynomial spectral trace to each F322W2 integration and extracted fluxes in each wavelength bin using an extraction aperture with a half-height of seven pixels centered on the estimated trace in each detector column. In doing so, we accounted for partial pixel effects in both the spectral and spatial directions.

tshirt

We reduced the data using the Time Series Helper & Integration Reduction Tool tshirt[15] (https://github.com/eas342/tshirt), which uses modified steps of the jwst pipeline to calculate rate images and performs a custom co-variance weighted optimal extraction\citeAppschlawin2020jwstNoiseFloorI.We followed a very similar procedure as previous NIRCam reductions[11]^,\citeAppahrer2022nircam.We begin by running the initial data quality, saturation and superbias subtraction steps with default parameters and reference of the jwst package\citeAppjwst_v1.10.2 version 1.10.2 and CRDS version 11.16.22.For Observation 8 (using the F322W2 filter), CRDS context jwst_1137.pmap, which we re-ran 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 the row-by-row odd/even by amplifier (ROEBA) step\citeAppschlawin2023defocused 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 setting the jump step’s rejection threshold to 6$\sigma$.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 co-variance weighted extraction\citeAppschlawin2020jwstNoiseFloorI 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 to 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 to a theoretical error (from photon and read noise) of 971 ppm.We note that the F444W has better 1/$f$ subtraction from sky pixels on 2 amplifiers whereas the F322W2 only has sky pixels available one one amplifier. The noise is elevated from theoretical expectations because there are still 1/$f$ correlations in the fast-read direction of the detector.

Eureka!

In Eureka!’s Stage 5, we fitted astrophysical and systematic noise models to our spectroscopic lightcurves. 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\citeAppstarry transit model with quadratic limb-darkening parameters fixed to an ATLAS model limb-darkening spectra\citeAppkurucz1993atlas9 as computed by the Exoplanet Characterization Toolkit (ExoCTK)\citeAppexoctk2021. We fixed all of our orbital parameters to those of Murphy et al.(in prep.; see the Derivation of Orbital Parameters section above) and adopted a minimally informative prior on the planet-to-star radius ratio.

Following the recommendations of ref.\citeAppbell2023_ers, we removed the first 800 integrations of the MIRI/LRS observations to remove the worst part of the initial ramp in the lightcurve. We also removed integrations 475–490 from the NIRCam/F444W observations which 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-Turns Sampler\citeApppymc3 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 the Gelman-Rubin statistic\citeAppGelmanRubin1992 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, while F444W exhibited only white noise only about 7% above the estimated photon limit on average. This excess noise is likely the result of insufficiently corrected 1/$f$ noise\citeAppschlawin2020jwstNoiseFloorI which gets partially converted to white noise in the time-series when spectroscopically binning the lightcurves. For MIRI/LRS, we used the estimated gain of 3.1 electrons/DN from ref.\citeAppbell2023_ers,bell2023arxiv, 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 $\mu$m but climbs steadily to approximately 50% above the limit near 10 $\mu$m and then climbs steeply to around 140% above the limit near 12 $\mu$m. Part of this excess long-wavelength white noise is likely due to the impact of background noise as these limits only considered the stellar photon noise. While neither of our NIRCam observations show evidence for residual red-noise in their Allan variance plots\citeAppAllan1966, several of our MIRI/LRS channels do show moderate amounts of residual red-noise. Following the $\beta$ error inflation method developed by ref.\citeAppWinn2008, we computed $\beta$, the ratio of the median value of our Allan variance plots at 15–23 minute timescales (around WASP-107b’s approximately 19-minute transit ingress/egress duration) 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 $\lesssim$6.5 $\mu$m and 0–30% for longer wavelengths. After this, we arrived at our final 2.45–12 $\mu$m transmission spectra.

Pegasus

We fit the F322W2 and F444W spectroscopic lightcurves from our Pegasus reduction using a BATMAN\citeAppbatman 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 prep.; 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\citeAppkurucz1993atlas9 by ExoCTK\citeAppexoctk2021 in the F322W2 and F444W bandpasses. To generate the limb-darkening coefficients we used the spectroscopic stellar properties measured in previous work[5].

To perform the lightcurve 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 maximization. We used the emcee\citeAppForemanmackey2013 Python package to perform the MCMC sampling, using twenty 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 the MCMC sampling had converged.

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

The estimated depth uncertainties were 1.15$\times$ the photon noise expectation in the F322W2 data and 1.06$\times$ the photon noise expectation in the F444W data.

tshirt

We followed a similar procedure as previous NIRCam light curve fitting[11]^,\citeAppahrer2022nircam.We use the starry\citeAppstarry code to model the lightcurves and pymc3\citeApppymc3 to sample the posterior distributions of the fits.We adopted uninformative priors on the stellar limb darkening with a 2 parameter quadratic law\citeAppkipping2013limbdarkening2param.This allows for comparison of results with the ATLAS stellar limb darkening models.The astrophysical starry model of the lightcurve is multiplied by an 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 lightcurve 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 No U-Turns sampling\citeApppymc3. We used 3,000 tuning steps, 3,000 sampling steps, and 2 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 lightcurve 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 panel a of Extended Data Fig.2). The mean depth uncertainty across the entire NIRCam wavelength range is approximately 90 ppm at our constant wavelength $\Delta\lambda=0.015$ $\mu$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}/\textrm{dof}=0.74$ for Eureka! vs. Pegasus, and $\chi^{2}/\textrm{dof}=0.64$ for Eureka! vs. tshirt. The larger difference between Eureka! and tshirt is largely because tshirt’s F322W2 spectrum consistently falls below that of Eureka! and Pegasus at wavelengths longer than $\sim$2.893 $\mu$m; this difference may be caused by tshirt’s ROEBA method since the deviation occurs right at an amplifier boundary (see panel a of Extended Data Fig.2). 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. Since the choice of reduction methods did not significantly impact 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.\citeAppkurucz1993atlas9,exoctk2021), (2) STAGGER-predicted quadratic law coefficients computed using ExoTiC-LD (refs.\citeAppmagic2015stagger,exoticld2022), (3) STAGGER-predicted 4-parameter law coefficients computed using ExoTiC-LD (refs.\citeAppmagic2015stagger,exoticld2022), and (4) freely fitting quadratic law coefficients using the reparameterization of ref.\citeAppkipping2013limbdarkening2param for efficient and physical sampling. As shown in panel b of Extended Data Figure 2, none of these spectra significantly differ from each other. In particular, we find that comparing the transmission spectrum of our fiducial analysis which 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 4-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 less 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 twenty-nine G102 and G141 spectroscopic transit lightcurves using a BATMAN\citeAppbatman transit model and standard detrending techniques for WFC3 timeseries data\citeAppbeatty2017kepler13. As with the JWST spectroscopic data, for these fits, we fixed the transit center time, orbital period, orbital inclination, the scaled semi-major axis, the orbital eccentricity, and the argument of periastron to previously measured values (Murphy et al., in prep.; 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\citeAppkurucz1993atlas9 by the Exoplanet Characterization Toolkit\citeAppexoctk2021. This left the planet-to-star radius ratio as the only free astrophysical parameter. We fit each spectroscopic channel individually.

All the G102 and the G141 lightcurves showed the usual ‘fishing-hook’ systematic trend within each orbit and the background linear slope typical for WFC3 timeseries observations. We modeled 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

Here, $t_{O}$ is the time since the start of each orbit’s observations, $t_{V}$ is the time since the start of the visit, and $A$ and $\tau$, and $m$ and $n$ are fitting coefficients for the exponential ramp and linear trend respectively. In addition to 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 lightcurve in each spectroscopic channel, we first performed a Nelder-Mead likelihood maximization to find an initial bestfit. We then used the emcee\citeAppForemanmackey2013 Python package to perform an MCMC exploration of the surrounding likelihood space to improve this bestfit estimate. In each channel we used twenty 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 below 1.1. We also checked to ensure that the median per-point flux uncertainty matched the standard deviation of the residuals to the bestfit lightcurve model in each channel. We additionally performed an Anderson-Darling test on each channel’s residuals to check for non-Gaussianity. The lightcurve residuals in each channel appear 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[2]. Instead, the updated 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 earlier analyses’ residuals. As a result, we do use the data from the third G141 orbit in our fitting. This slightly improves our measured depth uncertainties compared to those previous results.

MIRI

Unfortunately, while our 390 Hz noise removal and GLBS steps remove structured noise from the Lomb-Scargle periodogram\citeAppLomb1976,Scargle1982 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 TSO commissioning target L168-9b. For example, Extended Data Figure 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 $\mu$m bin sizes, but at that coarse of a resolution there are only 7 bins and there is substantial uncertainty in the estimated standard deviation from the small sample size. In addition, with finer bin size sampling than was used by ref.\citeAppbell2023_ers, we are better able to understand how the excess noise varies with varying bin size. In particular, while we find that the excess noise decreases with increasing bin size (similarly to what was reported by ref.\citeAppbell2023_ers), there also appears to be a spike in the excess noise around a bin width of 0.25 $\mu$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, while this noise spike seems to share an approximate period with the known 390 Hz noise, the 390 Hz noise removal step does not appear 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 likely 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 which 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 due to the minimal impact and high computational cost. That said, it does not appear that our 390 Hz noise removal procedure negatively impacts 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) to the European Consortium’s (hereafter referred to as the ‘EC’) three reductions presented by 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 by ref.[16], the transit depth differences among the EC spectra displays an increasing trend with wavelength with the CASCADe pipeline giving larger depths at longer wavelengths, which they attribute to different to different systematics models when fitting the lightcurves. As shown in Extended Data Figure 4, we find fairly similar conclusions with our fiducial spectrum which generally agrees well with the CASCADe spectrum. Looking more closely, it appears our fiducial spectrum has a slightly smaller transit depth around 7.5 $\mu$m and a slightly larger transit depths around 7 $\mu$m and from around 8–10 $\mu$m compared to the three reductions of ref.[16]. This results in our fiducial spectrum having a slightly smaller bump near the $\nu_{3}$ SO₂ feature and a slightly larger bump near the $\nu_{1}$ SO₂ feature.

1D-RCPE Grid-Retrieval

In order to produce physically self-consistent solutions we fit the transmission spectrum with a suite of models that impose 1D-radiative convective-photochemical-equilibrium (1D-RCPE). The coupling between the radiative-convective equilibrium solver (ScCHIMERA[11]^,\citeAppPiskorz2018, Mansfield2021) and the kinetics code (VULCAN[21]^,\citeAppTsai2017 using the H-O-C-N-S reaction list from ref.[21]) as well as additional details has been previously described in ref.[11]. The coupled model requires as inputs the incident stellar flux spectrum (T_eff=4,425 K, log($g$)[c.g.s.]=4.63; ref.\citeAppHusser2013) at the top of the planetary atmosphere (including the UV for photochemistry, using the HD 85512 MUSCLES spectrum\citeAppMusclesI,MusclesII,MusclesIII as a proxy, which matches the effective temperature and gravity of WASP-107 well within 1%) which can be scaled via an effective irradiation temperature, an internal temperature to set the deep adiabat, the atmospheric elemental abundances (via 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 Figure 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$\leq$[M/H]$\leq$2 in steps of 0.5), C/O (0.1$\leq$C/O$\leq$0.7 in steps of 0.2), internal temperature (100$\leq$T_int$\leq$500 in steps of 100 K), and a vertically constant eddy mixing (7$\leq$log${}_{10}K_{zz}$$\leq$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). Within 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, via 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 0.5$\leq$[M/H]$\leq$1.875 in steps of 0.125), C/O (0.05$\leq$C/O$\leq$0.6 in steps of 0.05), internal temperature (200$\leq$T_int$\leq$550 in steps of 50 K), and eddy mixing (7$\leq$log${}_{10}K_{zz}$$\leq$10 in steps of 0.5), resulting in 6,048 converged 1D-RCPE models. Based upon the temperature offset parameter from the course grid fit, this fine grid was generated at a cooler irradiation temperature (560 K) in order 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 1$\sigma$ (-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.

Within the Bayesian inference procedure, as described in ref.[11], we also fit for, on the-fly, the planetary reference pressure and reference radius (ref.\citeAppWelbanks2019a which shift and stretch the planetary spectrum and also affect the zero-point for the planetary gravity with height), an offset for the MIRI spectrum relative to the NIRCam F444 spectrum, the temperature profile offset parameter (described above), an ammonia abundance enhancement factor, a power law haze + gray 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 also now including SO₂, ref.\citeAppSO2_linelist) and fit to the data within 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]^,\citeAppBenneke2013. 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 significant contribution to the spectrum, as 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 Figure 5. We perturb the nominal 1D-RCPE solution (given in Extended Data Table 1) by 1) turning off mixing and photochemistry (equilibrium vs. 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 a $\sim$3 order-of-magnitude depletion in the methane abundance compared to 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₁₀CH₄ $<$ -5.6, details below). It is also apparent that T_int values lower than $\sim$350 K and a log${}_{10}K_{zz}$ lower than $\sim$8 struggle to result in the necessary reduction in the methane abundance even with mixing.

Cloud Parameterization

There are two leading modeling paradigms for modeling 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\citeAppGao+21_review. The later 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 as to unbias any inferences of other properties of interest (for example, volume mixing ratios of different gases or atmospheric metallicity, refs.[18]${}^{\text{,}}$\citeAppLine2016a,Welbanks2022,Barstow2020a). Some parametric models have attempted to capture the functional dependence of cloud extinction on wavelength using analytic models\citeAppTsiaras2018,Fisher2018 or Mie-theory\citeAppBenneke2019a, Pinhas2018, Welbanks2021, Grant2023.

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\citeAppAckerman2001, Molliere2019 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_{\rm 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 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_{\rm 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 SiO₂ best explained the MIRI spectral shape. However, the retrieved location of the cloud base was at relatively low pressures ($\sim$mbars)–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.\citeAppWakeford2015). While we are able to find satisfactory fits to the MIRI + HST observations (as in ref.[16]), when applied to the entire broad band spectrum with our NIRCam observations, we were unable to find satisfactory fits.

We also tested the plausibility of silicate clouds (MgSiO₃, SiO₂, and Mg₂SiO₄) and salt-sulfide clouds (Na₂S and KCl) in explaining the broad-band spectrum within the same eddysed\citeAppAckerman2001, Mai2019 framework. We use the 1D-RCPE atmospheric structure from an initial representative fit (T-P profile and chemistry using the nominal composition [M/H]=1.0, C/O=0.3, eddy diffusion coefficient–log$K_{zz}$=8.5, and a T_int=500 K) to self-consistently determine which clouds and where they form along the temperature-profile. We varied $f_{\rm sed}$ between 0.1 (vertically extended clouds) and 1 (compact clouds). Extended Data Figure 6 summarizes the condensate mixing ratios (panel a) and their impact on the spectrum (panel b) for an $f_{\rm sed}$=0.2. Only the lower $f_{\rm 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 of pressures. However, the lower $f_{\rm sed}$ values result in smaller droplets of 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 in-plausibility of the silicate cloud bases to exist at the altitudes probed by the observations that instead prompts us to take a phenomenological approach to the cloud modeling.

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 due to hazes as a scaling to Rayleigh-scattering[18]. 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 form

where $\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 functions; $\lambda_{0}$ is the wavelength at which the Gaussian is centered (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[18] and optically thick gray clouds as explained in refs.[11], in a linear combination with a cloud-free model to account for the presence of cloud/haze inhom*ogeneities\citeAppLine2016a. 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 employing 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 to 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.\citeAppWelbanks2022 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\citeAppWelbanks2021 and CHIMERA\citeAppLine2013b,Line2016a^,[2]. We select the later as our fiducial free retrieval comparison as 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 the ones 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 employed, 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 ${\sim}10^{-9}$bar to ${\sim}100$ bar, divided in 100 layers uniformly spaced in logarithmic pressure space. The vertical temperature structure is parameterized following the prescription from ref.\citeAppMadhusudhan2009. 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 assumes uniform mixing ratios for H₂O, CH₄, NH₃, CO, CO₂, SO₂, and H₂S using independent free parameters for each gas’ volume mixing ratio. Our choice of chemical species is motivated by those expected in exoplanetary atmospheres at these warm temperatures[9, 17]. As part of our initial analysis we considered the presence of PH₃ since 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 utilizing the one-sector parameterization of ref.\citeAppWelbanks2021 introducing the combined spectroscopic effect of a optically thick cloud deck with a cloud opacity $\kappa_{\text{cloud}}$, and hazes following an enhancement to Rayleigh-scattering[18] in a linear combination with a cloud-free atmosphere\citeAppLine2016a. 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\citeAppSkilling2004, Feroz2009, via PyMultiNest[45]${}^{\text{,}}$\citeAppFeroz2009,Feroz2019, 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: 7 molecular gases, 6 for the pressure-temperature structure of the planet, 1 for the reference pressure, 1 for the reference radius, 8 for the presence of inhom*ogeneous clouds and hazes, and 1 for an offset between the JWST MIRI observations and the JWST NIRCam F444W observations.

Our second free-retrieval with Aurora\citeAppWelbanks2021 follows largely the same model setup as the one explained above with CHIMERA. A detailed description of Aurora’s transmission modelling approach is available in ref.\citeAppWelbanks2021. 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.\citeAppRothman2010,Yurchenko2014,Yurchenko2011,Underwood2016; 2) the use of a free parameter P_cloud to account for optically thick clouds instead of $\kappa_{\rm{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 the CHIMERA retrieval at 68$\%$ confidence, and with predictions from the 1D RCPE grid retrieval. The results from this retrieval are also included in Table 2. Supplementary Information Figure 2 shows the equivalent to Figure 4 with the Aurora constraints, and highlights that while the inferences are consistent, the use of lower spectral resolution results in wider constraints in this case.

Supplementary Information Figure 1 shows the retrieved transmission spectrum and the contributions of the detected gases in WASP-107b’s atmosphere 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 in Figure 4 and Supplementary Information Figure 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 as shown in Extended Data Figure 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\sigma$ or greater. The reported detection significances in the main text result from the model comparisons with CHIMERA, while the detections with Aurora remain equally significant at $5\sigma$ 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\sigma$. Including wavelength dependent condensates is preferred at $13\sigma$ over just including optically thick clouds and Rayleigh-enhancement hazes. We note that the term ‘detection significance’ refer to a model preference between a reference model considering all gases and nested models for which each individual species is removed in turn\citeAppBenneke2013, Welbanks2021. 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\sigma$) that are difficult to interpret (see ref.\citeAppWelbanks2023 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 vs. 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.\citeAppVehtari2012, Vehtari2017, Vehtari2024) 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 out performs the other\citeAppWelbanks2023, McGill2023, Challener2023 and consequently, in this case, which wavelengths and instruments drive the NH₃ detection. Extended Data Figure 8 shows the difference in elpd ($\Delta$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$\sigma$ 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 the inclusion of NH₃ improved the model performance at all wavelengths (for example, there are points with positive $\Delta{\rm elpd}$ at all wavelengths), with a localized region of high performance data points near $3\,\mu$m. Indeed, the points with the highest preference for the model with NH₃, those with $\Delta{\rm elpd}>1$, are all between $2.9\,\mu$m and $3.2\,\mu$m and part of the NIRCam F322W2 observations. After NIRCAM F322W2, the next instrument with with positive $\Delta{\rm elpd}$ scores is HST-WFC3 G141 (9^th highest scoring point), followed by HST-WFC3 G102 (14^th highest scoring point) and NIRCAM F444W (seventeenth highest scoring point). MIRI’s highest scoring point is at 10.17 $\mu$m and corresponds to the 25^th highest scoring point with a $\Delta{\rm elpd}\sim 0.4$.

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

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

The NIRCam F322W2 provide the information necessary to confirm the tentative detection (${\sim}2$–3$\sigma$) 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 ${\sim}3.0\mu$m, and highlight an important advantage of NIRCam over other instruments on JWST that do not have the wavelength coverage or throughput necessary. While NH₃ has a strong spectral feature at ${\sim}10\,\mu$m, and the gas’ cross-section is dominant at $>8.6\,\mu$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 catalog magnitudes for the WASP-107 system to fit a model SED to the stellar emission. In conjunction with the Gaia\citeAppGaiaDR3 parallax for the system, this allowed us to estimate a stellar radius for the star WASP-107.

For our SED fits, we used catalog 2MASS JHK\citeApp2mass and AllWISE W1 and W2\citeAppwise 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 system’s parallax. We assumed $\log(g)=4.5$ and $[\mathrm{Fe}/\mathrm{H}]=0.0$ for the star WASP-107, which are both within 0.1 dex of the surface gravity and metallicity measured previously.[5] In practice, the precise values of $\log(g)$ and $[\mathrm{Fe}/\mathrm{H}]$ do not significantly affect the results from the SED fit.\citeAppstevens2018

We imposed a Gaussian prior on T${}_{\text{eff.}}$ based on previous spectroscopic measurements[5] of T${}_{\text{eff.}}=4425\pm 70$ K with the associated $1\sigma$ uncertainties as the prior width. Similarly, we imposed a Gaussian prior on the parallax to the system using the Gaia DR3\citeAppGaiaDR3 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\citeAppsf2011dustmap 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\citeAppbhac15 for the stellar SED. We computed a grid of surface luminosity magnitudes, corresponding to the bandpasses of the catalog magnitudes, for a range of T${}_{\text{eff.}}$ values. Since the BHAC15 models step by 200 K in T${}_{\text{eff}}$, we used cubic spline interpolation to estimate model magnitudes in between the points provided by the model atmospheres. We then scaled the interpolated surface magnitudes for the star by $R_{*}/d$ – where $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_{\odot}$. We then used the mean value of $R_{\rm p}/R_{*}$ as measured in our joint broadband transit fits to estimate the planetary radius of WASP-107b to be $R_{\rm p}=0.939\pm 0.019\,R_{\rm J}$.

Data Availability

The NIRCam data used in this paper are associated with JWST GTO program 1185 (PI 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 one-year exclusive access period. The MIRI data used in this paper are from JWST GTO program 1280 (PI Lagage; observation 1) and will also be publicly available on MAST at the end of their proprietary period. Additional intermediate results from this work are archived on Zenodo at https://doi.org/10.5281/zenodo.10780448 (ref.\citeAppwelbanks2024_wasp107b_zenodo).

Code Availability

We used the following codes to reduce and fit the JWST data: STScI’s JWST Calibration pipeline\citeAppjwst_v1.10.2, Eureka![14], tshirt[15], starry\citeAppstarry, PyMC3\citeApppymc3, numpy\citeAppnumpy, astropy\citeAppastropy2013,astropy2018, scipy\citeAppscipy, and matplotlib\citeAppmatplotlib.