Early Release Science of the exoplanet WASP-39b with JWST NIRCam

Measuring the metallicity and carbon-to-oxygen (C/O) ratio in exoplanet atmospheres is a fundamental step towards constraining the dominant chemical processes at work and, if in equilibrium, revealing planet formation histories. Transmission spectroscopy (for example, refs. 1,2) provides the necessary means by constraining the abundances of oxygen- and carbon-bearing species; however, this requires broad wavelength coverage, moderate spectral resolution and high precision, which, together, are not achievable with previous observatories. Now that JWST has commenced science operations, we are able to observe exoplanets at previously uncharted wavelengths and spectral resolutions. Here we report time-series observations of the transiting exoplanet WASP-39b using JWST’s Near InfraRed Camera (NIRCam). The long-wavelength spectroscopic and short-wavelength photometric light curves span 2.0–4.0 micrometres, exhibit minimal systematics and reveal well defined molecular absorption features in the planet’s spectrum. Specifically, we detect gaseous water in the atmosphere and place an upper limit on the abundance of methane. The otherwise prominent carbon dioxide feature at 2.8 micrometres is largely masked by water. The best-fit chemical equilibrium models favour an atmospheric metallicity of 1–100-times solar (that is, an enrichment of elements heavier than helium relative to the Sun) and a substellar C/O ratio. The inferred high metallicity and low C/O ratio may indicate significant accretion of solid materials during planet formation (for example, refs. 3,4,) or disequilibrium processes in the upper atmosphere (for example, refs. 5,6).

atmospheric metallicity of 1-100× solar (i.e., an enrichment of elements heavier than helium relative to the Sun) and a sub-stellar carbon-to-oxygen (C/O) ratio.The inferred high metallicity and low C/O ratio may indicate significant accretion of solid materials during planet formation e.g., 3,4 or disequilibrium processes in the upper atmosphere e.g., 5,6 .
JWST has demonstrated the necessary precision and wavelength coverage to make bulk characterization of hot exoplanet atmospheres routine 7 .The JWST Director's Discretionary Early Release Science (ERS) program provides the scientific community with observations of typical targets quickly enough to inform planning for the telescope's second cycle of scheduled observations.The primary goals of the Transiting Exoplanet Community ERS program (ERS-1366, led by N. M. Batalha, J. Bean, and K. B. Stevenson) are to demonstrate instrument capabilities, quickly build community experience, and seed initial discovery in transiting exoplanetary science 8,9 .The Panchromatic Transmission program observed a single exoplanet, WASP-39b, in transmission using four different instrument modes.It included overlapping wavelength coverage to cross-compare and validate all three near-infrared (NIR) instruments for time-series observations.The observations presented here form one quarter of this program, demonstrating the capacity of the Near-Infrared Camera (NIRCam) for transiting exoplanet atmospheric characterisation.WASP-39b is a highly inflated exoplanet of roughly Saturn mass, orbiting its G7 main-sequence star with a 4.05 day period 10 .We selected WASP-39b for its inactive host star and prominent spectroscopic features, which trace the atmospheric composition of the planet.We confirmed the star's relative inactivity through a photometric monitoring campaign using the Next-Generation Transit Survey (NGTS) 11 and Transiting Exoplanet Survey Satellite (TESS) 12 (see Methods).Reported atmospheric metallicities span a range of possible values (0.003 -300 solar) [13][14][15][16][17][18] due to limits on wavelength coverage, lower signal-to-noise ratio × data, and/or differences between analyses [19][20][21][22] .If the Solar System trend for gas giants 23,24 also applies to exoplanets, WASP-39b should have an atmospheric metallicity comparable to that of Saturn solar 25 ) and other Saturn-mass exoplanets.(10 × We observed a single transit of WASP-39b with JWST's NIRCam instrument on 22-23 July 2022 (19:28 -03:40 UT).The Grism R and F322W2 filter in the long-wavelength (LW) channel dispersed light from 2.420 -4.025 µm at a spectral resolution R of 1570 -2594 over 1023 resolution elements.The short-wavelength (SW) channel allowed the simultaneous measurement of light, i.e. photometry, spanning 2.0 -2.2 µm using the WLP8 weak lens and F210M filter.See the Methods section for more details.
The team conducted three independent reductions of the NIRCam LW spectroscopic data and four independent fits and analyses of the reduced data.We also performed two independent analyses of the SW photometric data.For both data reductions (LW and SW), customising the JWST Science Calibration Pipeline (jwst) to allow for minor adaptations to default steps and values worked best (see Methods).The wavelength solution available with the reference files provided by the JWST Calibration Reference Data System (CRDS) at the time of our analysis was inaccurate (particularly for the blue edge of the LW channel), so we redefined our wavelength values using a polynomial wavelength calibration derived from a planetary nebula observed as part of commissioning (Program 1076).The subtle variation in transit depth around 2.8 µm is due primarily to water vapour in the planet's atmosphere.The vertical striping in the residuals is due to 1/ noise.
We found no large systematic structures affecting the LW light curves and a minuscule ramp at the start of the SW light curve (Figure 1, panel d).The only other systematic identified was 1/ noise (or pink noise, where is frequency), which describes the detector's correlated read   noise 27 .For NIRCam, this manifests as weak structures in the dispersion direction, as shown in Figure 1c.We did not correct for 1/ noise in the final LW reduction because it did not  impact the precision reached by individual spectroscopic light curves (compare tshirt and Eureka! in Figure 2 for analyses with and without 1/ noise corrections).We removed  structures due to 1/ noise in the SW reduction (see Methods).We found a linear model in  time was sufficient to detrend the data, which produced uncertainties 1.18 the photon noise × limit (median of 135 ppm for the transit depths) at a binned spectral resolution of 15 nm ( 15 ∼ pixels).Similarly, the photometric transit-depth precision was 1.35 the noise limit at × 53 ppm.The residuals are Gaussian (see Extended Data Figure 5 in Methods).

Figure 2:
The transit spectrum of WASP-39b as measured from JWST's NIRCam instrument.The coloured points with 1 uncertainties depict our independent analyses of the σ spectroscopic LW channel (2.420 -4.025 µm) and photometric SW channel (2.0 -2.2 µm) with their respective throughputs shown in grey.All analyses agree with the broadband Spitzer point (black circle, 3.2 -4.0 µm).The broad feature centred at 2.8 µm spans 2.5 scale heights ( 2,000 km) and is due primarily to water vapour within WASP-39b's atmosphere.
∼ We note the consistency between analyses in the fine structure.
Figure 2 displays the independently derived transit spectra and photometry.Each reduction is consistent with our selected reduction (Eureka!) to better than 1 , as is the broadband 3.6 σ µm Spitzer point 13 .The overall shape of the spectrum is due primarily to absorption of water vapour (feature centred at 2.8 µm).The right-axis scale is in equivalent scale heights, where one scale height is approximately 800 km.
To interpret the presence of other molecules within the planetary atmosphere, we compared the Eureka!transit spectrum to a set of independently computed atmospheric model grids that spanned a range of cloud properties, metallicity values, and C/O ratios (see Methods).  ) 13 .We also see  3) leads to a more tentative identification here.Each 2 forward model grid prefers significant cloud coverage, which impacts the spectrum at mbar ∼ pressures, despite differing cloud parametrisations between grids with varying levels of physical complexity (see Methods).
In a hot ( 1000 K) solar-metallicity atmosphere with a stellar C/O ratio, CH would be ∼ 4 visible as a strong peak at 3.3 µm (grey dashed line in Figure 3, Extended Data Figure 7 in Methods) under thermochemical equilibrium.Such a peak is absent in the reduced spectrum.We quantified this using a residual fitting test (see Methods).In a higher-metallicity and/or lower-C/O atmosphere, carbon is increasingly partitioned into CO and CO , and the CH above which the goodness of fit per free parameter, , gets increasingly worse (i.e., 2).
We also tested whether other data reductions favoured best-fit models with stronger methane abundances, but found they did not have any statistical significance.Driven by this CH upper limit, the single best fit from each grid favours the lowest C/O ratio 4 (0.229, 0.3, and 0.35 for PICASO 3.0, PHOENIX, and ATMO, respectively) within that grid.These best-fit point values for C/O from the three grids agree well with the value of found by Ref. 13 .We examined the effect of an even lower C/O grid point by 0. 31 −0.05 +0.08 computing the best-fit PICASO 3.0 model with a C/O of 0.115, but found no discernible difference in the transit spectrum.Comparing our inferred C/O ratio for WASP-39b's atmosphere to that of its host star, we see that it is substellar ( 0.35, whereas WASP-39 is ≤ 0.46 0.09, 23 ).We also note that the C/O ratio shown here represents the carbon-to-oxygen ± fraction of the planet's upper atmosphere rather than that of the whole atmosphere, as these NIRCam observations probe approximately the 0.1 -10 mbar pressure range.WASP-39b's temperature-pressure profile is cool enough for the formation of silicate (i.e., O-bearing) cloud species at depth, which would deplete oxygen from the upper atmosphere and actually increase the C/O ratio aloft compared to the bulk planetary envelope 28,29 .Figure 4 compares our best-fit metallicity values, shown as separate O and C abundances, and C/O ratios to previous studies using HST data, as well as results for exoplanets observed at high resolution and Solar System gas giants.The JWST/NIRCam data rule out a super-stellar C/O ratio for WASP-39b.Additionally, Figure 4 demonstrates JWST's capability to measure the C/O ratios of giant planet atmospheres by observing both O-and C-bearing species , which until now has only been achieved through high-resolution exoplanet observations e.g., 30,31 .Similar measurements have been difficult to achieve from HST alone..Blue ≥2σ 2 points show all previous estimates of the metallicity of WASP-39b from HST data, offset in mass for clarity 13,[15][16][17][18] 23 .Our results for WASP-39b favour a super-stellar volatile abundance and substellar C/O ratio.However, we emphasize that a full retrieval will be necessary to determine accurate means and error bars for the NIRCam results.1σ The apparent substellar C/O ratio inferred from chemical equilibrium models may trace photochemical processes in the planet's upper atmosphere.For example, photochemical destruction of CH in the upper atmosphere could explain the absence of a CH peak at 3.3 µm e.g., 6,34 .The most likely immediate products of CH photolysis, such as HCN or C H , would be produced in abundances too small ( a few parts-per-million (ppm), 6,34 ) to be ≲ robustly detected with a single NIRCam transit, even from complete CH conversion.

4
Alternatively, much of the carbon available from CH photolysis could have been oxidized by 4 photodissociated H O to form CO and CO 6,[34][35][36] , though the absolute abundances of these two carbon reservoirs would not have been meaningfully altered since their abundances under chemical equilibrium are already higher than that of CH .Other proposed disequilibrium 4 chemistry processes could reduce the CH abundance at the terminator without also decreasing the C/O ratio 5,[37][38][39][40] .We defer the exploration of complex disequilibrium models to atmospheric retrieval analyses using the full set of data provided by the Transiting Exoplanet Community ERS program.That dataset will also constrain the presence of additional oxygenand carbon-bearing species to provide a more robust constraint on the C/O ratio than we can obtain here.However, the C/O ratio estimate we report from NIRCam is broadly consistent with the C/O ratio found from the other individual ERS WASP-39b datasets, which range from best-fits that are sub-solar (NIRISS/SOSS, Feinstein et al., submitted; NIRSpec/PRISM 3.0 -5.0 µm, 7 ; NIRSpec/G395, Alderson et al., submitted) to a slightly super-solar upper limit (NIRSpec/PRISM 0.5 -5.5 µm, Rustamkulov et al., submitted).
If disequilibrium chemistry is not prevalent in the planet's upper atmosphere, the inferred high metallicity and low C/O ratio can be tied back to WASP-39b's formation.The most prominent scenario is that WASP-39b formed via core accretion exterior to the water-ice line and accreted low-C/O solid material in situ and/or while migrating inward within the protoplanetary disk 4,45,46 .Taken as such, JWST observations could offer important clues regarding the degree to which hot Jupiter atmospheres undergo solid accretion during their early evolution.
Here, we have demonstrated the excellent performance of NIRCam for exoplanet transmission spectroscopy.With the first JWST exoplanet spectra now comparable to the first near-infrared Jupiter spectra 47 , the future promises many exciting discoveries and major advancements in the formation, evolution, and atmospheric chemistry of hot Jupiters.

Methods
As part of this article's Reproducible Research Compendium, located on Zenodo at https://doi.10.5281/zenodo.7101283,we provide a Jupyter Notebook with step-by-step data reduction instructions replicating our chosen analysis, saved outputs from various pipeline stages, and the data used to generate relevant figures.

Photometric Monitoring of Host Star
In order to confirm that WASP-39 is a relatively inactive star, and that the JWST observations were not adversely affected by stellar activity, we carried out photometric monitoring with the ground-based Next Generation Transit Survey (NGTS) 11 .Monitoring began at the end of April 2022 and continued until late August, spanning the JWST ERS transit observations of WASP-39b in July.We used one camera on most photometric nights to take a series of 10 s images lasting on average for 2 h.The resulting monitoring light curve is plotted in Extended Data Figure 1 (top panel), showing one binned point for each night.Also included is the TESS sector 51 PDCSAP (Pre-search Data Conditioned Simple Aperture Photometry) light curve of WASP-39 12 , which is binned to 2 h to be comparable with NGTS.Both light curves have been detrended against sky brightness.They show evidence for stellar activity, but only with a low amplitude of 0.06% in NGTS.Also plotted in Extended Data Figure 1 (lower panels) are individual transit observations of WASP-39b with NGTS and TESS (the times of which are indicated on the monitoring light curve).For four of the NGTS transits, we employed multiple cameras.This significantly improves the photometric precision 48 , which is otherwise limited by atmospheric scintillation 49 .The transit models were generated from the system parameters listed in Extended Data Table 1.We only fit the transit times and the mutual depth of the TESS transits, which is slightly shallower than expected.

Extended Data
The transit observations in Extended Data Figure 1 show no evidence for starspot-crossing events, which would be visible as bumps in the transit light curve.The absence of such events across multiple high-precision transits provides additional evidence that WASP-39 is a quiet star and that the JWST ERS transit observations are unlikely to be adversely affected by stellar variability.

JWST NIRCam Observation
JWST observed the 2.8-hour transit of WASP-39b over a span of 8.2 hours, providing a baseline before and after transit to measure transit depths accurately.A dichroic beam splitter allows NIRCam to simultaneously observe a target in both short wavelength (SW) and long wavelength (LW) channels 50,51 .The LW channel used the Grism R + F322W2 filter to observe a wavelength range of 2.420 -4.025 µm with a spectroscopic resolving power of at 4 µm (Extended Data Figure 2, top panel).The SW imaging channel used the ∼1, 600 WLP8 weak lens and F210M filter (2.0 -2.2 µm) to produce the hexagonal pattern shown in the bottom panel of Extended Data Figure 2. Spreading the light prevents saturation, reduces variability due to image motion over an imperfect flat field, and allows monitoring of mirror-segment alignment.Both SW and LW channels used the SUBGRISM256 subarray mode with four output amplifiers and the SHALLOW4 readout pattern to minimize data volume.With 12 groups per integration (82.17 s total), we acquired 366 integrations for this transit observation.No impactful tilt events were noted in this observation.

Data Reduction and Calibration
We conducted independent data analyses using multiple pipelines and fitting tools to ensure that we obtained the same transmission spectrum using different reduction pipelines.We also varied the fitting methods within a given data reduction pipeline.
Many of the reductions presented below used intermediate data products from or made minor edits to the JWST Science Calibration Pipeline (jwst) 1 , which we briefly summarize here.jwst is a Python software suite for processing data from all JWST instruments and observing modes, and is organised into three stages.Stage 1 takes in uncal.fitsfiles and performs detector-level corrections and ramp fitting for individual exposures (i.e., ramps-to-slopes conversion; these ramps are the flux increases during an exposure, not to be confused with baseline ramps over the course of the entire transit).Stage 2 takes in slope images (ramps) from Stage 1 and performs an assignment of the world coordinate system, flat fielding, and assignment of a wavelength solution.Stage 3 takes in calibrated 2D images from Stage 2 and extracts a time series of 1D spectra.The default pipeline settings include a flux calibration step at Stage 2. In all data reductions presented below, we skipped that step, as it introduced scatter in the extracted spectral time series.This is justified because the transit depths we compute are relative, rather than absolute, flux measurements.
Below we describe the independent data reductions applied to the SW photometry and LW spectroscopy, respectively.In each case we note where data reductions deviated from the standard jwst pipeline.

Short-Wavelength Photometry
We performed two independent short wavelength data reductions using the open-source Eureka! and tshirt pipelines.

Eureka! Short Wavelength Reduction
Eureka! is an open-source pipeline designed to perform spectral extraction and fitting for JWST exoplanet time-series observations 52 .The Eureka! short wavelength data reduction used the default jwst settings for Stages 1 and 2, with the exception of increasing the rejection threshold during jump detection to , which improved the quality of the 10σ resulting light curve.
In Stage 3, we first masked all pixels for which the "DO_NOT_USE" data quality flag was raised by the jwst pipeline.We then performed an outlier rejection along the time axis for each individual pixel in a segment using a threshold, repeating this process twice.Next, 7σ we corrected for the 1/ noise in each of the four amplifier regions by subtracting the median  flux in each row calculated without pixels containing the star.We interpolated over flagged pixels using a cubic function.Finally, we determined the image centre and performed aperture photometry on the target.We explored different target apertures and background annuli, and chose the combination that minimised the root-mean-square variations, leading to a target aperture radius of 65 pixels and a background annulus from 70 to 90 pixels relative to the centre.

tshirt Short-Wavelength Reduction
tshirt is an open source pipeline 2 that has tools to modify jwst and performs photometric and optimal spectral extraction of light curves.
In the Stage 1 SW analysis, tshirt applied a row-by-row, odd/even-by-amplifier (ROEBA) subtraction algorithm that used background pixels to reduce the 1/ noise.In this  procedure, background pixels are used to correct each group in a similar fashion as with reference pixel correction. 3The ROEBA correction happens after the bias subtraction step.First, the median of all even columns' background rates is subtracted from all even columns and the median of all odd columns' background rates is subtracted from all odd columns to remove most pre-amp reset offsets and odd/even pixel effects.Next, the median of each column's background rate is subtracted from each row to remove the 1/ noise for timescales  longer than a row read time (5.24 ms).The correction was applied to each group so that 1/ noise would not be detected as spurious jumps or cosmic rays by the pipeline.We used all pixels more than 201 pixels from the source to estimate the background and 1/ noise, then  subtracted the median of each row from all pixels in that row.Stage 2 of jwst was skipped, as it only changes the rates from ADU per second to physical units and conducts flat fielding.This does not affect the relative measurements of the light curve (due to the high pointing precision) and allows for comparison to detector-level effects.
For the photometric extraction, we used a source radius of 79 pixels and a background annulus of 79 to 100 pixels.We performed a 2D Gaussian fit to determine the centre of the aperture.

Long-Wavelength Spectroscopy
We performed three independent long-wavelength data reductions, using the Eureka!, HANSOLO, and tshirt pipelines.
The reference files in the Calibration Reference Data System (CRDS) at the time of our analysis included a linear solution for wavelength as a function of coordinate (the  dispersion direction), but this is not strictly accurate at the blue end.For all methods, we use commissioning program 1076 to derive a third-degree polynomial wavelength solution that uses the Pfund and Bracket Hydrogen Series in the planetary nebula IRAS 05248-7007.The residuals in this solution are nm and the stellar absorption lines in WASP-39 agree with ≲0. 1 the solution to within 1 nm.The difference between the corrected wavelengths and the original wavelength solution is almost zero at the red end of the spectrum, but increases to about 50 nm at the blue end.

Eureka! Long-Wavelength Reduction
We investigated several variations of the Eureka!long-wavelength data reduction to minimise the MAD of the final extracted light curves, with different settings for cosmic-ray jump detection, identifying the spectral trace, the aperture size for spectral extraction, the region for background subtraction, and limits for outlier rejection.Here we present details of the data reduction that produced the spectrum shown in the main body of the paper.
Stages 1 and 2 were identical to the jwst pipeline, with the exception of increasing the rejection threshold during jump detection to .In Stage 3, we first trimmed the data to a 6σ subarray extending from pixels in the cross-dispersion direction and in 4 − 64 4 − 1704 the spectral direction.We then masked any pixels with NaN values for the flux or error.We fit the spectral trace with a Gaussian profile and corrected for the curvature of the trace to the nearest integer pixel.We excluded a region 14 pixels wide on either side of the spectral trace from the background calculation and performed a column-by-column linear fit to subtract the background.We used a double-iteration threshold for outlier rejection of the sky 7σ background along the time axis during background subtraction.Additionally, we used a 7σ threshold for outlier rejection during the polynomial fit to the background.To obtain the spectrum, we constructed a normalised spatial profile using the median of all data frames, then used optimal extraction 53 on an aperture with a half-width of 9 pixels.For the optimal extraction, we rejected outliers above a threshold.Extended Data Figure 3 shows the 10σ curvature-corrected, background-subtracted, median frame with indicated background and aperture regions.
Extended Data Figure 3: Curvature-corrected, background-subtracted, median frame.We perform optimal spectral extraction on the pixels in between the green dashed lines.We use the pixels outside of the two orange solid lines for background subtraction.The flux spans -200 -1000 electrons, thus drawing attention to the residual background features.Right: Vertical slice depicting the flux averaged over detector pixels 855 to 865.The background region clearly demonstrates some low-level residual structure.

HANSOLO Long-Wavelength Reduction
The HANSOLO (atmospHeric trANsmission SpectrOscopy anaLysis cOde) pipeline was originally developed to analyse ground-based transmission spectra observed with 8m-class telescopes 54,55 and was adapted to enable its use on NIRCam data.HANSOLO begins with the calibrated rateints.fitsoutputs of jwst Stage 1.
We used the LACOSMIC algorithm 56 to remove cosmic ray effects from the two-dimensional images and identified the spectral trace using a Moffat function fit to each column.To remove the sky, we fitted and subtracted a linear trend from each column, excluding from the fit a region of 20 pixels on either side of the trace centre.We then extracted the spectrum by summing over an aperture with a half-width of 3 pixels.The spectra from different images were aligned with each other using cross correlation.To correct outlier pixels, each spectrum was normalised to account for the effect of the transit on the time series.Outliers away from the mean were removed from the time series of each > 3σ wavelength point in the normalised spectra and replaced with the median value over time.We then rescaled the spectra to their original amplitudes.

tshirt Long-Wavelength Reduction
As with the short-wavelength reduction, a few modifications were made to the Stage 1 CalWebb ramps-to-slopes pipeline.ROEBA subtraction reduced 1/ noise (described above  for photometry); however, only pixels 1847 to 2044, which are on the rightmost amplifier, are available as low-illumination background.
For Stage 3, tshirt performed optimal spectral extraction weighted by the covariance between pixels 27 .We used a spectral aperture centred at pixel 34 in the spatial direction with a half-width of 5 pixels.We selected the background region to extend between pixels 5-24 and 44-65 in the spatial direction.The background was fit with a column-by-column linear trend with 3-sigma clipping.For the spectral extraction, we fit the spatial profile with a cubic spline with 20 knots and an outlier rejection threshold of 30 .If a σ pixel was deemed an outlier either by the "DO_NOT_USE" data quality flag or by the spatial profile outlier detection, the rest of the spatial profile was weighted by the reference profile to ensure that the flux was conserved.For the covariance weighting, a correlation of 8% was assumed between pixels as measured by background pixels' noise properties.

Data Analysis and Fitting
We used both Eureka! and tshirt to fit the short-wavelength light curves.In both cases, the light curves were fit with models that included both the transit and systematic noise.However, in order to investigate the effect of different systematic models on the resulting spectra, each fit used a slightly different noise model.Extended Data Table 1 summarises the systematics models which were used in each short-wavelength fit.
For the long-wavelength fits, we summed the data into 15 nm bins ( pixels).We ∼ 15 experimented with bins as small as 10 nm, but found that reducing the bin size below 15 nm led to poor constraints on the limb darkening and added additional scatter to the resulting spectrum.Extended Data Figure 4 shows that the noise is primarily Gaussian out to long time scales of order the length of ingress/egress.Additionally, we created a white light curve by summing the extracted spectra over the entire 2.420 -4.025 µm wavelength region.We experimented with different wavelength cut-offs but chose to extract spectra in this wavelength region because the low instrument throughput affected the quality of the extracted light curves beyond this region.Extended Data Figure 5 shows all reduced transmission spectra with one bin added on the blue end and two added on the red end, as well as the relative throughput at the wavelengths of these bins.This figure demonstrates the large error bars derived from data near the edges of the NIRCam/F322W2 bandpass.Therefore, we recommend that future works limit extracted spectra to the wavelength region between 2.420 -4.025 µm.1: Details of the two methods used to fit the short-wavelength photometry.Abbreviations for priors are as follows: U=uniform prior, with numbers indicating lower and upper limits; N=normal prior, with numbers indicating mean and sigma; LN=log-normal prior, with numbers indicating mean and sigma.

Extended Data Table
Extended Data Figure 4: Normalised root mean square error as a function of bin size for all spectroscopic channels.The red line shows the expected relationship for perfect Gaussian white noise.The black lines show the observed noise from each spectroscopic channel for the Eureka!long wavelength reduction.Values for all channels are normalized by dividing by the value for a bin size of 1 in order to compare bins with different noise levels.The black lines closely follow the red line out to large bin sizes of ( -hr time ≈30 ≈0. 5 scales), which demonstrates that the residuals to the fit are dominated by white Gaussian noise.
Extended Data Figure 5: Top panel: Transmission spectra from our reductions when including additional data on the blue and red edges (now spanning 2.405 -4.055 µm), demonstrating the large error bars and diverging data points near the edges of the NIRCam bandpass in the LW spectroscopic channel.The differences in retrieved transmission spectra by subtracting the Eureka!spectrum from the other three reduced spectra shown in the top panel.This shows the strong agreement between the spectra; however, we do note minor disagreements at shorter wavelengths that we attribute to differences in the treatment of limb-darkening effects within the individual fitting methods.
We fit the long-wavelength light curves using four independent pipelines: chromatic-fitting, Eureka!, HANSOLO, and tshirt.chromatic-fitting is an open-source 4 Python tool to perform light-curve fitting, built on the data visualiser chromatic (Berta-Thompson, in prep. 5).For this work, chromatic-fitting light-curve fitting was applied to a Eureka!data reduction.As with the short-wavelength fits, we fit the long-wavelength light curves with models that include different noise parameterisations.Extended Data Table 2 summarises the systematics models that were used in each long-wavelength fit.2: Details of the four fitting methods used to fit the long-wavelength spectroscopy.Abbreviations for priors are as follows: U = uniform prior, with numbers indicating lower and upper limits; LU = log-uniform prior, with numbers indicating lower and upper limits; N = normal prior, with numbers indicating mean and sigma; LN = log-normal prior, with numbers indicating mean and sigma.The notation ``spec-fixed'' indicates that a value was fit in the white light curve and fixed to the best-fit value for the spectroscopic light curves.* Note that different fitting methods used different parameterisations of the planetary orbit and noise model, so not all methods fit for all of the listed parameters.Parameters marked with ``N/A'' were not fit in that method and were instead derived from the other fitted parameters.† Subscripts 0, 1, etc. indicate the 0th, 1st, etc. order terms in polynomial models.‡ The GP model was only applied to the white light curve.For the spectroscopic light curves, the divide-white method62 was used to remove the GP systematics.For all fits, the parameters were estimated with a Markov Chain Monte Carlo (MCMC) fit, using either the emcee Python package 58 (for fits performed with Eureka!), the pymc3 Python package 59 (implemented through the Exoplanet code 60,61 , for fits performed with chromatic-fitting or tshirt), or the CONAN Python package 54,55 (for fits performed with HANSOLO).The number of free parameters and the resulting differential MADs of the light curves from each fit are also listed in Extended Data Tables 1 and 2. The best-fit parameters from the white light-curve fits are given in Extended Data Table 3.

Extended Data Table 3:
Best-fit orbital parameters from both short-wavelength (SW) and white-light long-wavelength (LW) fits.
In the process of performing the fits to the long-wavelength data, we regularly found that the best-fit transmission spectra were shifted vertically for different limb-darkening parameterisations and, for some reductions, exhibited changes in the apparent size of the water feature.In particular, we found that light-curve fits with all limb-darkening coefficients fixed to outputs from ExoTiC-LD [63][64][65] could result in a biased planet spectrum and might present a higher level of time-correlated noise in the residuals.We attribute this to a combination of JWST's high-precision light curves and deficiencies in the stellar limb-darkening models to accurately represent WASP-39 66,67 .Therefore, the results presented here use the quadratic limb darkening law, in its classical form or reparameterised by Ref. 57 , with one or both coefficients as free parameters.We confirmed that these parameterisations produce transmission spectra that are consistent both with each other and with the spectra resulting from using more complex limb darkening parameterisations, such as a four-parameter law with either fixed or free parameters 68 .We therefore recommend that future transmission spectrum analyses with NIRCam use similar methods.Limb-darkening conclusions from the full Transiting Exoplanet Community ERS program will be discussed further by Espinoza et al. (in prep.).The final fitted light curves are shown in Extended Data Figure 6 and the final transmission spectra are shown in Figure 2.Both the short wavelength and long wavelength datasets are also available in our Reproducible Research Compendium on Zenodo at https://doi.10.5281/zenodo.7101283.The median difference between each transmission spectrum and the Eureka!spectrum is (using the maximum error at each point), 0. 87σ which demonstrates a remarkable level of agreement.Additionally, the residuals showed no evidence for time-correlated noise, as shown in Extended Data Figure 5.
For ease of interpretation, we compared our atmospheric models to only one transmission spectrum.We selected the Eureka!spectrum, as it was on average nearest the median spectrum (the median transit depth at each bin).
Each model was binned to the resolution of the data to perform the fitting.We performed χ 2 these three independent model grid fits to fully vet our inferences about the atmospheric metallicity and the presence of specific molecular features within the data.
Extended Data Figure 6: Top Row: Time-series NIRCam data for the WASP-39b system, from three independent spectral extractions.Color represents relative brightness at each time and wavelength, normalized by the median stellar spectrum.Bottom Row: Resulting residuals after fitting the time-series NIRCam data.
The PICASO 3.0, Vulcan, & Virga Model Grid Our primary atmospheric model grid is built from the open-source radiative-convective equilibrium code PICASO (Planetary Intensity Code for Atmospheric Spectroscopy Observations) 75 , version 3.0 76 , which was developed from the Fortran-based EGP [77][78][79] .We used PICASO 3.0 to generate one-dimensional temperature-pressure profiles in thermochemical equilibrium.The base PICASO 3.0 forward model grid computes atmospheric mixing ratios using variations of planetary intrinsic temperatures (T ) of 100 K,  200 K, and 300 K; C/O ratios of 0.229, 0.458, 0.687, 0.916; and atmospheric solar metallicity values of 0.1, 0.316, 1.0, 3.162, 10.0, 31.623,50.119, and 100 solar.The PICASO grid × assumes full day-night energy redistribution.To compute model transmission spectra from the atmospheric profiles, we used opacities described by Ref. 79 (see in particular Extended Data Table 2); which sources H O from Refs. 80,813][84] , CO from Ref. 85 , and H S from Refs. 82,86,87. 2 We then used the one-dimensional CHON-based chemical kinetics code VULCAN 34 and the cloud modelling code Virga 88 , which is the Python implementation of the Eddysed cloud code 89 , to post-process disequilibrium chemistry from mixing and photochemical products as well as the effect of clouds.These additional post-processed grids also include vertically constant eddy diffusivities (K ) of 10 -10 cm /s in steps of 2 dex, and both clear and cloudy models.For the Vulcan disequilibrium runs, we only computed model grid points for a select subset of metallicity values (1, 10, 50, and 100 solar) and C/O ratios × (0.229, 0.458, 0.687).We found that neither the cloudy nor clear disequilibrium grids from VULCAN offered an improvement in the value.Given the sparseness of these pre-computed χ ν 2 disequilibrium grid models, we left rigorous quantification of self-consistent disequilibrium chemistry in the atmosphere of WASP-39b to future work.
Within PICASO, clouds are implemented both as grey absorbers and as Mie scatterers using temperature-relevant cloud condensate species from Virga.For the grey clouds, the grid specified a cloud optical depth ( ) between 1 and 0.1 bar ranging from = 3.2 τ  τ  × 10 to 1 in steps of 0.1 dex across all wavelengths.For clouds of specific condensates, we

The PHOENIX Model Grid
We also used a grid of atmosphere models from the PHOENIX radiative-convective equilibrium code to fit the data [72][73][74] .Similar to the PICASO 3.0 grid, parameters including the day-night energy redistribution factors, interior temperature (200 K, 400 K), bulk atmospheric metallicity (0.1, 1, 10, 100 solar), and C/O ratio (136 grid points from 0.3 to 1) × were varied.Aerosol properties were parameterized through a haze factor (0, 10 multi-gas × Rayleigh scattering) and a grey cloud deck pressure level (0.3, 3, and 10 mbar).Models with molecular abundances quenched at 1 bar to simulate vertical mixing were also calculated.The grid also included rainout to account for species sequestered as condensates in the deep atmosphere.Opacities are described by Refs. 74,90and taken from Ref. 86 .

Grid Fits to JWST/NIRCam Data
We applied each of our three grids -ATMO, PHOENIX, and PICASO 3.0 -to fitting the NIRCam F322W2 spectrum (2.4 -4.0 µm).In doing so, we found that the models strongly favoured a solar-or super-solar-metallicity atmosphere (1 -100 solar), a sub-stellar C/O × ratio ( 0.35), and substantial contribution from clouds, which are parameterized differently ≤ by each model grid (see each grid description above).We show the best fits from each model grid in Extended Data Figure 7.This interpretation is in agreement with the results using JWST's NIRSpec/PRISM instrument from 3.0 -5.0 µm 7 , improving on the wider spread from previous HST-only [13][14][15]17,18,92 or HST and ground-based optical interpretations 16 . The bestfit equilibrium model from the PHOENIX grid had 100 solar metallicity, a C/O × ratio of 0.3, and a cloud deck at 3 mbar.Cloudy models were generally preferred over clear models, but not with statistical significance ( of 1.25 compared to 1.22). Th PHOENIX grid χ ν finds best fits with very high metallicity (100 solar), so this low confidence regarding × clouds reflects the cloud-metallicity degeneracy inherent in data restricted to narrow wavelengths e.g., 93 , as well as potentially the sparseness of the model grid.

Extended Data
For the ATMO grid, the best-fit equilibrium model to the NIRCam spectrum was 1 solar × metallicity, a C/O ratio of 0.35, a cloud factor of 5 and a haze factor of 1.As with the other two grids, strongly cloudy models (cloud factor of 5) were preferred to clear models ( of ≥ χ In Extended Data Figure 8, we show the comparison between the spectra of HST/WFC3 (G141 and G102, covering 0.8 -1.65 µm) and JWST/NIRCam (F210M+F322W2, 2.0 -4.0 µm).We chose to only show Wide Field Camera 3 (WFC3) observations from HST, as these are of higher precision than observations from the Space Telescope Imaging Spectrograph (STIS) or ground-based data 13 .Additionally, as HST/WFC3 has the most archival exoplanet data of any instrument on Hubble, future JWST exoplanet programs will primarily rely on this HST instrument for inter-telescope comparisons or extending the wavelength coverage of JWST data.For example, the addition of optical and shorter wavelength near-infrared data can help break metallicity degeneracies by better constraining the presence and extent of clouds 13,e.g., 93 .High altitude clouds or hazes can be inferred from their particle sizes, where small particles scatter shorter wavelengths more efficiently e.g., 94,95 , thus enabling the disentanglement of a very cloudy, low metallicity atmosphere from a less cloudy, high metallicity atmosphere 17 .

Molecular Detections
Once we found the "single best fit" for the PICASO grid to the NIRCam spectrum (10 solar, × C/O = 0.229, grey cloud optical depth = 2.6 10 from 1 to 0.1 bar), we used this as a base × −3 model to explore the significance of specific molecular detections.First, we tested whether we could improve the best fit in the presence or absence of H O, CO , CH , or H S. We reran the best-fit base model by zeroing out each of these species in turn, shown in Figure 3, and then repeating our analysis.χ 2 We found that while the presence of H O, H S, and CH resulted in a better value, only H O and H S did so in a statistically meaningful way.Since H S does not contain strong molecular features within the NIRCam wavelength range, the Gaussian residual fitting we perform for the detection significance of other molecules is not applicable, and we left its further quantification to more rigorous atmospheric retrieval analyses.Increasing the CH show the residual features left after subtracting out the gas in question (CO , top, and H O, With the best fit in hand, we investigated the presence of individual molecular species.For molecular detection significances, we performed the same Gaussian residual fitting, shown in Extended Data Figure 9, as for the detection of CO in the NIRSpec/PRISM 3.0 -5.0 µm 2 analysis 7 .We find a Bayes factor, ln( ), of 123.2 between the Gaussian residual and constant  models for H O over the whole NIRCam wavelength range, corresponding to 15.9 , a strong 2 σ detection.For CO we find ln( ) of 0.82 between the Gaussian residual and constant models 2  between 2.4 and 2.9 µm, or 1.9 , which is a weak or non-detection 96 .CO is strongly σ 2 detected at 4.3 µm in the NIRSpec data for WASP-39b Gaussian residual fitting for CH and find a weak or non-detection at approximately 2 .the atmosphere of WASP-39b must rely on wavelengths that can fully capture the complete absorption feature, which is beyond the reach of high fidelity NIRCam/F322W2 measurements.

Figure 1 :
Figure 1: The relative brightness of the WASP-39 planetary system as a function of time and wavelength, as measured by NIRCam.Panels a, b, and c show the spectroscopic data while panels d, e, and f present the photometric (short wavelength) channel.The extracted flux normalised by the median stellar spectrum is shown in panels a and d, the best-fit transit and systematic models are shown in panels b and e, and the residuals are shown in panels c and f.The flux decrease results from the transit of exoplanet WASP-39b in front of its star.The subtle variation in transit depth around 2.8 µm is due primarily to water vapour in the planet's atmosphere.The vertical striping in the residuals is due to 1/ noise.

Figure 3
Figure3shows a representative best-fit model highlighting the contributions of major molecular absorbers.

Figure 3 : 4 ( 2 4 2 σ
Figure 3: Contributions of key absorbers impacting the spectrum.Top: The best-fit PICASO 3.0 equilibrium model (10 solar, C/O = 0.229, moderate grey clouds with cloud × optical depth of 2.5 10 ) is shown compared to the Eureka!reduction, along with models × −3 with individual molecular species removed to show its contribution to the spectrum.Each model is normalized to the data for illustration by offsetting each model to have the same transit depth at 2.8 µm.Water predominately sets the shape of the spectrum, followed by the influence of clouds.The grey dashed line shows a cloudy solar-metallicity and stellar-C/O atmospheric model, illustrating the lack of a strong CH peak seen in the data.Bottom: The 4 opacities of the dominant molecular species at an optical depth of 1 in the atmosphere.In the single best-fit model shown in the lower panel, the methane peak at 3.3 µm is blended out by water absorption.However, manual scaling of CH gives an upper limit of CH abundance 4 4(blue line) for the single best-fit model shown in the top panel.

2 4 peak at 3 . 3
µm disappears.Therefore, the absence of a strong CH peak at 3.3 µm in our data 4 drives the metallicity to higher values and the C/O ratio to lower values.We scaled the CH 4 volume mixing ratio (VMR) within our single best-fit PICASO 3.0 model (10 solar × metallicity; C/O ratio of 0.229) to determine an upper limit on the abundance of CH at 1 4 mbar, where it contributes most strongly to the spectrum.Within our single best-fit model scaling, we find an upper limit on CH abundance at 1 mbar of 5.5 10 (or 55 ppm) VMR,

Figure 4 :
Figure 4: Trends in elemental abundances and C/O ratio with planet mass.Panels a, b, and c show the abundances of oxygen, carbon, and net volatiles (oxygen + carbon), respectively, scaled to stellar values.Grey points in panel (a) show HST constraints based on H O detections, with the grey dashed line showing the best-fit trend from Ref.18 .Blue ≥2σ

Figure 1 :
Photometric monitoring of WASP-39 (top) and individual transit observations (bottom) using NGTS (magenta) and TESS (dark purple).The black marks indicate the times of the four JWST ERS transit observations.The monitoring light curve shows evidence for optical variability, but with an RMS amplitude of only 0.06% in NGTS.The times of the individual transit observations are indicated on the top panel, and they are all consistent with transits free of starspot crossings or other features associated with stellar activity.

Figure 2 :
Raw NIRCam image of the LW (top) and SW (bottom) channels.The faint horizontal stripes seen in the LW channel originate from neighbouring objects.The SW channel is able to track changes in alignment for individual mirror segments.

− 6 used
Virga to compute log-normal particle size distributions using sedimentation efficiency (f ) values of 0.6 to 10 for MnS, Na S, and MgSiO along the range of K .Smaller efficiencies, f , with larger eddy diffusivities, K , generated more extended   cloud decks and stronger cloud opacity.

Figure 7 : 2 1. 22 .χ ν 2 resulted 9 . 3 ≥
Measured transmission spectrum compared to atmospheric forward model grids.Top: The single best fit for each model grid (shown as solid colored lines; PICASO 3.0, ATMO, PHOENIX), fits the planet spectrum (Eureka!reduction) with χ ν All single best fits prefer at least solar metallicity and substantial cloud cover.Also ≤ shown as a grey dashed line is a solar metallicity, stellar C/O ratio atmospheric model, demonstrating the lack of methane absorption seen in the spectrum.Because we can put an upper limit on the CH abundance, the preferred C/O ratio found by the model grids is 4 sub-stellar.Bottom: Residuals of each best fit, shown as the model spectrum subtracted from the reduced spectrum and divided by the uncertainty in transit depth.The residuals show wavelength-dependent correlations, the origin of which are unknown and left for a future study.For the NIRCam-only fit, the PICASO grey-cloud scheme produced a slightly better best fit ( = 1.16) than the PICASO + Virga more realistic clouds ( = 1.23), both of which were χ clear-model best fit (100 solar) with = 1.53.The Virga best-fit grid × in an atmosphere of 1 solar metallicity, C/O = 0.229, f = 0.6, and K = 10 cm/s ×   This Virga best-fit model consists of clouds of MnS and MgSiO with deep ( 100 bars) 2 cloud bases and diminishing optical depth up to mbar pressures.∼

8 :
Our JWST/NIRCam spectrum compared to existing HST/WFC3 data.As in Extended Data Figure 7, but with the addition of HST/WFC3 data from 0.8 to 1.65 m, showing the comparable precision and complementary wavelength µ coverage offered by the combination of NIRCam and HST/WFC3.

4Extended Data Figure 9 :
abundance beyond that of the best-fit model also improved the , though again not to high χ Gaussian residual fitting of H O and CO .The blue points 2 2

2 2 bottom
) from the single best-fit model.The Gaussian model ensemble fit to the residual is shown in red; the best-fit Gaussian ensemble to a flat-line model is shown in blue.We strongly detect H O at nearly 16 and show weak evidence for CO (small feature at 2.6 µm)

4 σ 2 points ( 4 .
Both WASP-39b NIRSpec datasets7, Rustamkulov et al., submitted, Alderson et al., submitted  observed evidence for a molecular feature near 4.0 µm, which is currently best explained by SO .The reddest data 025 µm) from NIRCam also show an increase that is consistent with this feature > seen in the NIRSpec data.However, as shown in Extended Data Figure5, these NIRCam data points have very large error bars because the detector throughput drops off dramatically past 4.0 µm.Future investigations to thoroughly explore the physicochemical likelihood of SO in 2 NIRSpec/PRISM at 4.3 µm 7 , but the overlap between the CO feature at 2.8 µm and the 2 − 1. 37 −0.13 +0.05 weak evidence for CO absorption, previously seen with high confidence using 2 2 broad H O feature (illustrated in Figure
7, Rustamkulov et al., submitted, Alderson et al., submitted, but the strong overlapping H O band at 2.8 µm prevents NIRCam from making a significant 2 CO detection.Given our upper limit on CH abundance, we also performed the same