Identification of carbon dioxide in an exoplanet atmosphere

Carbon dioxide (CO2) is a key chemical species that is found in a wide range of planetary atmospheres. In the context of exoplanets, CO2 is an indicator of the metal enrichment (that is, elements heavier than helium, also called ‘metallicity’)1–3, and thus the formation processes of the primary atmospheres of hot gas giants4–6. It is also one of the most promising species to detect in the secondary atmospheres of terrestrial exoplanets7–9. Previous photometric measurements of transiting planets with the Spitzer Space Telescope have given hints of the presence of CO2, but have not yielded definitive detections owing to the lack of unambiguous spectroscopic identification10–12. Here we present the detection of CO2 in the atmosphere of the gas giant exoplanet WASP-39b from transmission spectroscopy observations obtained with JWST as part of the Early Release Science programme13,14. The data used in this study span 3.0–5.5 micrometres in wavelength and show a prominent CO2 absorption feature at 4.3 micrometres (26-sigma significance). The overall spectrum is well matched by one-dimensional, ten-times solar metallicity models that assume radiative–convective–thermochemical equilibrium and have moderate cloud opacity. These models predict that the atmosphere should have water, carbon monoxide and hydrogen sulfide in addition to CO2, but little methane. Furthermore, we also tentatively detect a small absorption feature near 4.0 micrometres that is not reproduced by these models.

WASP-39b is a hot (planetary equilibrium temperature of 1170 K assuming zero albedo and full heat redistribution), transiting exoplanet that orbits a G7-type star with a period of 4.055 days 15 .The planet has approximately the same mass as Saturn (M = 0.28 MJup) but is ~50% larger (R = 1.28 RJup), likely due to the high level of irradiation it receives from its host star [16][17][18] .We chose this planet for the JWST ERS transmission spectroscopy observations because analyses of existing space-and ground-based data detected large spectral features and showed that there was minimal contamination of the planetary signal from stellar activity 10,[19][20][21] .The main spectral features previously detected were confidently attributed to sodium, potassium, and water vapor absorption 10,[19][20] , while CO2 was suggested to explain the deep transit at 4.5 µm seen with Spitzer 10 .Atmospheric metallicity has long been thought to be a diagnostic of the relative accretion of solids and gas during the formation of gas giant planets, both of which bring heavy elements to the hydrogen-dominated envelope and visible atmosphere [4][5][6] .The metallicity of WASP-39b's host star, which is a proxy for the metal enrichment of the protoplanetary disk that the planet formed in, is approximately solar 15,[22][23][24] .Therefore, the planet mass -atmospheric metallicity trend observed in the solar system giants 25,26 predicts that it has an enhancement of ~10x solar (like that of Saturn, Ref. [27]).Additionally, interior structure models that match WASP-39b's low density predict a 95th percentile upper limit for the atmospheric metallicity of 55x solar, under the limiting assumption that the planet has no heavy element core and that all the metals are evenly distributed throughout the envelope 28 .
Despite having some of the highest signal-to-noise detections of spectral features in its transmission spectrum, modeling of the existing data for WASP-39b has resulted in metallicity estimates ranging across five orders of magnitude, from 0.003x to 300x solar 10,[29][30][31][32][33] .The wide range of values stems from the data being of insufficient quality to break the degeneracy between clouds and metallicity in transmission spectra models 34 , as well as uncertainty over the interpretation of the photometric measurements by the Spitzer Space Telescope at 3.6 and 4.5 µm.Thus, spectroscopic data with greater precision, finer spectral channels, and wider wavelength coverage were needed to better constrain the metallicity of this (and other) giant exoplanet atmospheres. .The first JWST ERS observation of WASP-39b was obtained using the Near Infrared Spectrograph (NIRSpec) 35,36 on July 10, 2022, between 15:24 and 23:37 UTC.We used the Bright Object Time Series (BOTS) mode with the 1.6" x 1.6" fixed slit aperture and the PRISM disperser to capture spectra between 0.5 and 5.5 µm.The data were recorded using the SUB512 subarray with five groups per integration and the NRSRAPID readout pattern, which gave integration times of 1.38 s.NIRSpec obtained a total of 21,500 integrations over 8.23 hours of observations centered on the 2.8 hour transit duration of WASP-39b.
The count rate in the PRISM mode varies significantly over the bandpass due to the spectral energy distribution of the star and the wavelength-dependency of the spectrograph dispersion.Therefore, the observations were designed to saturate at shorter wavelengths in order to obtain sufficient signal-to-noise ratio at the longer wavelengths in the bandpass for first-of-its kind spectroscopy of transiting exoplanets.Wavelengths between 0.71 and 2.09 µm have at least one group saturated in the pixel at the center of the spectral trace.We concentrate here on the analysis of the data longward of 3.0 µm that are not impacted by saturation to investigate the spectrum overlapping with the previous 3.6 µm and 4.5 µm Spitzer photometric measurements.The subset of the PRISM data described herein has a native spectral resolving power (R=/Δ, where  is wavelength) of 100 -350.For this study, we binned the data to lower resolving powers (values range from 60 to 200 depending on wavelength and reduction).The binning is done at the light curve level before the fitting of the transit depths that constitute the transmission spectrum.Analyses of JWST/NIRSpec transit observations obtained during commissioning have shown that similar levels of binning as we use here results in minimal systematics 37 .An analysis of the complete PRISM dataset at full resolution including recovery of the saturated part of the spectrum is ongoing.We reduced the NIRSpec PRISM data for WASP-39b using the JWST Science Calibration Pipeline along with customized routines to minimize noise in the time series spectra (see Methods).We performed four different reductions of the transmission spectrum starting from the uncalibrated data 21,[38][39][40] .Figure 1 shows derived spectroscopic transit light curves from one of the reductions.We confirm with our analysis of the WASP-39b data that NIRSpec transit observations at a resolving power of 60 -200 are nearly free of systematics.We achieved close to photon-noiselimited measurements in the spectroscopic light curves after trimming the first 10 minutes of data and removing a linear trend in time with an average rate of ~190 ppm/hour across the bandpass.We also obtained similar results by fitting the full time-series with a downward trending exponential ramp (timescale ~100 minutes) combined with a quadratic function of time.The lack of large systematics in these data stands in contrast to previous transit spectroscopy observations with space-or ground-based telescopes 41 .
Figure 2: Independent reductions of the WASP-39b transmission spectrum.The JWST data (small colored points) are compared to Spitzer's two broadband photometric measurements (grey circles and corresponding sensitivity curves labeled "IRAC1" and "IRAC2").The axis on the right shows equivalent scale heights (750 -1000 km) in WASP-39b's atmosphere; for plotting purposes we assume that one scale height corresponds to 800 km.The JWST data are consistent with the Spitzer points (within 2σ) when integrated over the broad bandpasses (indicated by the horizontal lines).The relative transit depths between the 3.6 and 4.5 µm channels are also consistent within 2σ between independent reductions of the JWST data, with most of the deviation coming from the 3.6 µm bandpass.Vertical error bars indicate 1σ uncertainties.
The transmission spectra derived from the different reductions, shown in Figure 2, have excellent agreement.They all exhibit a large feature at 4.3 µm, as well as a smaller feature near 4.0 µm (discussed below).Detailed modeling of the FIREFLy-reduced data yields a statistical significance scale he ghts of 26σ for the large feature (see the Methods).We attribute this feature to CO2 absorption based on a comparison of the resolved band shape to theoretical models and the spectra of brown dwarfs 42 .Figure 2 also includes Spitzer's two broadband photometric measurements 10 , which are consistent with the JWST data to better than 2σ after integrating the transmission spectrum over the Spitzer bandpasses.We also see good agreement (better than 2σ for all reductions) in the relative transit depths between the 3.6 and 4.5 µm channels.The comparison shown in Figure 2 demonstrates both the consistency in the derived spectra from multiple, independent analyses and the reliability of the previous Spitzer measurements.The key parameters of the model are 10x solar metallicity, carbon-to-oxygen ratio of 0.35, and cloud opacity of 7x10 -3 cm 2 /g.The impact of the opacity sources expected from thermochemical equilibrium over the full bandpass are indicated by removing the opacity contribution from individual gasses one at a time.As in Figure 2, the axis on the right shows equivalent scale heights in WASP-39b's atmosphere.The bottom panel shows the molecular absorption cross-sections for each gas in the best-fit model.The model is well matched to the data ( !/ "#$# =1.3), suggesting that our assumptions broadly capture the important physics and chemistry in WASP-39b's atmosphere.However, there is a feature near 4.0 µm that cannot be reproduced by the models used here.The strong CO2 absorption (4.1 -4.6 µm) and the apparent lack of methane (3.0 -3.5 µm) is what drives the solution to an elevated atmospheric metal enrichment, ruling out previous low metallicity estimates [29][30][31] .The other reductions and models give similar results.We compared the data to a suite of one-dimensional atmospheric structure and transmission spectrum models to constrain the composition of WASP-39b's atmosphere.These models assume radiative-convective-thermochemical equilibrium, and they adopt a scaled solar abundance pattern.We calculated planet-specific grids of these models over a range of atmospheric metallicities, carbon-to-oxygen ratios, and cloud properties using four different codes.These grids of self-consistent model transmission spectra were then fitted to the FIREFLy-reduced data (the fit results are independent of which data set we use) while also adjusting for a reference radius at 1 bar.The results are illustrated in Figure 3; see the Methods for further details.
Under similar assumptions, all four model grids are able to match the dominant spectral morphologies -namely the strong CO2 feature between 4.1 and 4.6 µm and the rise in transit depth blueward of 3.6 µm due to water vapor (a species that had been detected previously at shorter wavelengths, Ref. [10]).More subtle modulations over the whole bandpass are potentially due to contributions from clouds, carbon monoxide, and hydrogen sulfide, though the degree to which the two gas species contribute is unknown pending further study.
Several models for warm gas giant atmospheres predict that the CO2 abundance scales quadratically with atmospheric metallicity, becoming detectable at 4.3 µm for metallicities above that of the Sun 1-3 .The representative best-fit model shown in Figure 3 is consistent with this scenario.It has a 10x solar metal enrichment and a slightly sub-solar carbon-to-oxygen ratio (0.35, compared to the solar value of 0.55; Ref. [43]).The moderate contribution of cloud opacity predicted by the best-fit model is consistent with interpretations of previous population-level studies of planets that have similar temperatures and gravities as WASP-39b 44,45 .It is also consistent with the predictions of aerosol microphysics and global circulation models of hot giant planets [46][47][48] .
In addition to the large CO2 feature, we also identify a smaller spectral feature near 4.0 µm that is not matched by our thermochemical equilibrium models (see Figure 3).This feature is present in all four independent reductions and has a significance of 2σ (see the Methods).Further data analysis and modeling including non-equilibrium chemistry are needed to fully assess the robustness of this feature and to identify the chemical species that gives rise to it.Additional JWST ERS observations of WASP-39b that will use the G395H grating on NIRSpec also have the potential to confirm the 4.0 µm feature and resolve it in greater detail.
The grid fits explored here favor lower metallicities than Refs.[10, 21], and higher metallicities than Ref [31], even though the Spitzer data that their studies included are consistent with our JWST data.The higher precision and more resolved measurement of the CO2 feature enabled by JWST pulls the models of Refs.[10, 21] to lower metallicity and increased cloudiness.Nevertheless, it is not possible to obtain a robust confidence interval on this inference without more rigorous Bayesian analyses, which is left to future work (see Methods).Continued modeling of WASP-39b will also be aided by the future measurements of the planet's transmission spectrum from 0.5 -5.5 µm that are also being obtained by this Early Release Science program.The final transmission spectrum will ultimately have higher spectral resolution than the data presented here (>4x over most of the bandpass) and will be validated using multiple JWST instruments.

Data reduction
We reduced the JWST NIRSpec PRISM data for WASP-39b using four separate pipelines to confirm that the results did not depend on the specifics of the analyses, as was sometimes the case for results from the Spitzer Space Telescope (e.g., Ref [49]).The descriptions below refer to calibration pipelines and other software whose code and citations appear in the Code Availability section, below.

tshirt pipeline
We used the Time Series Helper and Integration Reduction Tool 40 (tshirt) to extract light curves of the spectrum.This pipeline modifies the JWST Calibration pipeline steps to improve the precision of the reduction.tshirt has been used to successfully analyze the JWST transit observations of HAT-P-14b that were obtained during commissioning with NIRCam 37 .First, we used an updated bias frame from commissioning program 1130 observation 29 and ran the JWST Calibration pipeline until the reference pixels step.We then applied a correction for 1/f-noise which varies for odd and even rows and for each column.We use background pixels for the calibration since reference pixels are not available in this subarray.We skipped the jump and dark subtraction steps because they were seen to add noise to the light curves.tshirt fits the profile of the spectrum with splines and rejects outlier pixels that are more than 50 sigma from the spline fits.We used covariance-weighted extraction 50 with an assumed pixel correlation of 0.08.For spectral extraction, we used a background region no closer than 7 pixels on either side of the source and an extraction region width of 16 pixels.The scatter in the light curve was consistent with the theoretical limit of photon and read noise over short timescales.
We fit the light curves with a second-order (quadratic) polynomial baseline, uninformative quadratic limb darkening priors, and an exponential startup ramp with 10σ clipping of outliers.To begin, we fit the white light curve with priors on the transit center, inclination, and period from Ref [22].We also used the a/R* from Ref [22] but widened the uncertainty on this parameter because the enforced prior resulted in significant residuals.Next, we fit each spectroscopic light curve individually with the orbital parameters fixed at the value from the white light posterior medians.We modeled the light curves using the "exoplanet" code 51 and the pymc3 52 sampler.We evaluated the wavelengths using the JWST Calibration pipeline at pixel row 16 (Y=16) from the world coordinate solution.This uses an instrument model and could not be verified due to a lack of strong stellar absorption features at the NIRSpec resolution.All the other reductions adopted this wavelength calibration.As shown in Figure 1, the standard deviation in the out-of-transit light curve approaches the theoretical limit of photon and read noise at short wavelengths, but is 20% to 40% higher at longer wavelengths, which may be related to uncorrected 1/f-noise.

Eureka! pipeline
Eureka! 39 is a data reduction and analysis pipeline for time-series observations with JWST or HST.Its modular, multi-stage design provides flexibility and ease of comparison at any step, starting from uncalibrated FITS files and resulting in precise transmission or emission spectra.Eureka! has been used to successfully analyze the JWST transit observations of HAT-P-14b that were obtained during commissioning with NIRCam 37 .
We began the data reduction process using the UNCAL files available from the MAST archive.The first stage of the Eureka!pipeline is primarily a wrapper for Stage 1 of the JWST Calibration pipeline, which converts groups to slopes.For this dataset, we skipped the jump detection step as it led to a large fraction of detector pixels being incorrectly flagged as outliers.We did, however, search for and flag outliers at multiple points in subsequent stages.We also manually updated the bad-pixel map to include identified hot pixels on the detector that were not provided in the current (July 2022) full-detector STScI data quality map.As part of Eureka!, we performed a custom background subtraction at the group level prior to Stage 1 ramp fitting to account for 1/f-noise introduced during detector readout.We set the top and bottom six rows of the detector as our background region and flagged pixels deemed outliers at >3σ.We then subtracted the mean flux per pixel column and repeated this for each group and integration in the observation.Similar to Stage 1, the second stage of the Eureka!pipeline is a wrapper for Stage 2 of the JWST Calibration pipeline, which calibrates the 2D time series of fitted slopes.Here, we skipped the flux calibration step, thus leaving the data in units of DN/s.
For Stage 3, we performed background subtraction and optimal extraction of the stellar spectrum for each integration with Eureka!.We only used pixels 14 to 495 in the dispersion direction of the 512x32-pixel subarray, as NIRSpec's throughput is negligible beyond this range.We also masked pixels that have a non-zero data quality flag to avoid any impact of outlier pixels on the extracted spectra or background subtraction.The position of the source on the detector along the cross-dispersion dimension is located by fitting a Gaussian to the pixel values summed over all detector columns.For each pixel, we examined its flux variation in time and performed a double-iteration, 10σ outlier rejection test.We then executed a second column-by-column background subtraction, this time at the integration level, using pixels located at least 8 pixels away from the source position to compute the mean background per column.Performing this additional background subtraction reduced the number of outliers in the measured light curves and accounted for the residual background and/or noise introduced during the ramp fitting procedure.As with Stage 1, we exclude 3σ outliers from our background region.We adopted an aperture halfwidth of 7 pixels for our optimal spectral extraction step, constructing the profile from the median frame.At the end of this stage, we obtained a time series of 1D spectra.
For the remaining stages, we used multiple pipelines (Eureka! 39and ExoTEP [53][54][55] ) to generate and fit the light curves.We first generated median-normalized light curves at the instrument's native resolution (i.e., from each detector column) using our Stage 3 outputs.We then clipped additional outliers in time for the white and spectroscopic light curves.For this step, we first rejected integrations that were more than 3σ outliers for the source position in the crossdispersion direction, the width of the fitted Gaussian to the spatial profile, or the drift in the dispersion direction.Next, we produced a median-filtered version of the light curve and clipped out 3σ outliers in flux.We jointly fit astrophysical and systematics model parameters to the white and individual spectroscopic light curves.Our astrophysical transit model used the batman package 56 with uniform priors, fitting for the following astrophysical parameters: the two coefficients of a stellar quadratic limb-darkening law, impact parameter, semi-major axis, transit time, and the planet-to-star radius ratio (a/R*) in each of the wavelength channels.While the limbdarkening coefficients and planet-to-star radius ratio were fit independently in each spectroscopic channel, we used the best-fitting value of the planet's impact parameter, semi-major axis, and transit time from a white light curve fit as a fixed value in the wavelength-dependent fits.For the systematics model, we assumed a linear trend in time for each wavelength channel, fitting for both the slope and y-intercept.Last, we fit a single-point scatter to each light curve, which illustrates the level of additional noise required for our joint model to reach a reduced chi-squared of unity.The white light curve residuals have an RMS of 3013 ppm and the spectroscopic light curves above 3 µm have a median RMS of 5779 ppm.Similar to the reduction shown in Figure 1, both pipelines reach near photon noise.The Eureka! and ExoTEP transmission spectra appear nearly identical; therefore, only one (Eureka!) is shown in Figure 2.

Tiberius pipeline
We built upon the pipeline developed for the analysis of LRG-BEASTS data 21,57,58 to provide an independent reduction of the data.We began with the outputs of the JWST Calibration Stage 1 pipeline with the jump step correction turned off.We created bad-pixel and cosmic-ray masks by identifying 5σ outliers in running medians operating along pixel rows and along individual pixels in time.Prior to tracing the spectra, we interpolated each column of the detector onto a finer grid, 10x the initial spatial resolution, in order to improve the extraction of flux at the sub-pixel level.We used a 4th-order polynomial to trace the spectra and a 4-pixel-wide aperture.To remove the 1/f-noise, we fit a linear polynomial to 21 background pixels along each column in the cross-dispersion direction.Next, to correct for shifts in the dispersion direction, we crosscorrelated each stellar spectrum with the first spectrum of the observation to account for very small (0.003 -0.005) sub-pixel shifts.Our white light curve spans a wavelength range of 0.518-5.348μm after masking saturated pixels, and our 147 spectroscopic light curves used 3-pixel-wide bins across this same wavelength range.We masked frames 20751-20765 due to a high gain antenna move that led to increased noise in the light curves.
We fit our light curves with a combination of a quadratically limb-darkened transit model (through batman 56 ) with a linear-in-time polynomial.We began by fitting the white light curve to derive the system parameters: inclination, i, time of mid-transit, TC, the semi-major axis scaled to the stellar radius, a/R*, and the linear limb darkening coefficient, u1.We placed wide boundaries on the parameter values only to prevent unphysical values.In practice, the parameter values did not get close to the boundaries.We fixed the planet's orbital period to 4.0552941 days and eccentricity to 0 from Ref [22].We fixed the quadratic coefficient, u2, to theoretical values determined by Exo-TiC-LD 59,60 with 3D stellar models 61 , and fit for u1.We used a Levenberg-Marquardt algorithm to fit our light curves, rescaled our photometric uncertainties to give a reduced  2 =1 for our best-fit model and then re-ran the fits.For the spectroscopic light curves, the system parameters (i, TC, a/R*) were held fixed to the best-fit values found from the white light curve.The white light curve residuals had an RMS of 2761 ppm and the spectroscopic light curve residuals had a median RMS of 6731 ppm.In both cases, the variance of the residuals scales upon binning as expected for Poisson noise.

FIREFLy pipeline
We also reduced the data using the Fast InfraRed Exoplanet Fitting Lyghtcurve (FIREFLy) reduction routines 38 .These routines utilize the JWST Calibration pipeline with custom modifications.This pipeline has been used to successfully analyze the JWST transit observations of HAT-P-14b that were obtained during commissioning with NIRSpec G395 37 .We removed 1/fnoise (see Ref [36]) at the group level, as the 1/f-noise changes from group to group.We also skipped the jump step and instead flagged and removed cosmic rays, bad pixels, hot pixels, and other outliers using median filtering of the data both spatially and in time, flagging pixels using a 5σ outlier threshold algorithm.The time series of 2D spectra were aligned using cross-correlation and interpolation, with the time series spectra exhibiting an RMS jitter of 0.005 pixels in the Xaxis direction and 0.0026 pixels in the Y-axis direction.We found a small inverse ramp in the light curves, which settled down after the first 2000 exposures, which we discarded.We fit the light curves with the batman 56 transit model along with a linear baseline and a second-order jitter detrending polynomial of X and Y detector position as described by Ref [38], which are present in the spectrophotometry at the 53+/-2 ppm level in the x-direction and 140+/-3 ppm in the ydirection.We applied a fixed quadratic limb darkening law using the 3D models 61 computed using the methods of Ref [62] from ExoTiC-LD 59,60 .In fitting the 3 to 5.5 μm white light curve, we allowed the semi-major axis in units of stellar radii a/R*, inclination i, and central transit time T0 to freely vary along with the transit depth and systematics model.We used the Markov chain Monte Carlo sampling routine EMCEE 63 to find the best-fit parameters and measure the posterior distribution.We find the 3-5.5 μm white light curve has a transit depth of 2.1368+/-0.0014%and achieves 808 ppm scatter in the residuals.This is within 6% of the expected noise limit of 758 ppm as calculated by the JWST Calibration pipeline, with the scatter of the residuals decreasing to below 40 ppm upon binning with no detectable red noise.We fit each spectroscopic light curve shown in Figure 2 with the same astrophysical and systematic models as the white light curve, except fixing the system parameters (a/R*, i, T0).The transmission spectral light-curve residuals for each bin are typically within 5% of pipeline error or better, also with no detectable red noise.

Data-Model Comparison
We compared the extracted transmission spectral data to a suite of 1D self-consistent radiativeconvective-thermochemical equilibrium model atmospheres (see e.g., Refs [64, 65] for a general description of such models) described below.In short, all models are able to fit the 3-5.5 μm spectra consistently (with  2 / "#$# < 1.4) with a 10x solar metal enrichment and varying grey cloud opacity for their single best estimate.Comparisons of the model fits from each grid are shown in ED Figure 1.For additional parameters within the grid (e.g., C/O and heat redistribution), there is some discrepancy between each model grid's single best estimate values.Additional Bayesian analyses are needed to rigorously quantify confidence intervals on atmospheric properties of interest, which is beyond the scope of this work.Future works will focus on modeling that includes the effects of disequilibrium chemistry, aerosol microphysics, and three-dimensional circulation effects.We assumed the following parameters in the modeling: stellar Teff = 5512 K, stellar radius = 0.932 RSun, planet mass = 0.281 MJup, planet radius = 1.279RJup, and planet orbital semi-major axis = 0.04828 AU.
ED Figure 1: Comparison of transmission spectrum modeling results from different codes for WASP-39b.Despite different radiative-convective equilibrium and chemical solvers, treatments of clouds, grid spacing, and grid-fitting approaches, all four grids arrive at the same 10x solar metallicity point solution.Additionally, all four provide an acceptable fit to the data, with best fitting  !/ "#$# < 1.4.

ScCHIMERA
This framework was first described in Refs.[66, 67], with the most recent updates, methods, and opacity sources described in Ref. [68].We compute the converged atmospheric structure (temperature-pressure and thermochemical equilibrium gas mixing ratio profiles) over a grid of atmospheric metallicity ([M/H], where brackets indicate log10 enrichment relative to solar 43 ) spaced at 0.25 dex intervals between 0 and 2.25 (1 to 175 times solar) and carbon-tooxygen ratio (C/O) at values 0.2, 0.35, 0.55, 0.7, 0.75, 0.8.We assume full day-to-night temperature redistribution 69 as planets in this temperature regime are unlikely to possess strong day-to-night temperature contrast 70,71 .We then compute transmission spectra 72,73 from these converged atmospheric structures.To match the models to the data, the DYNESTY 74 fitting routine is used to search for the optimal [M/H] and C/O (via nearest neighbor) while also simultaneously adjusting the 1 bar planetary radius (which controls the absolute transit depth) and an opaque, grey, uniformly vertically distributed, cloud opacity ( %&" ).The optimal model resulting from this process is [M/H] = +1.0,C/O = 0.35, and log10cld = -2.15cm 2 /g.The metallicity and cloud opacity are primarily driven by the strength of the 4.3 µm CO2 feature and lack of CH4 absorption near 3.3 µm.This result is what is shown in the main text (see Figure 3), which also illustrates the relative contribution of the key opacity sources (H2O 75,76 , CO 77,78 , CO2 79,80 , H2S 78,81 , and CH4 78,82 ) to the overall spectral shape.ED Figure 2 shows the atmospheric structure (temperature profile and gas mixing ratio profiles) for this best fit model.2: Atmospheric structure arising from the best fit model.The thick red curve (and corresponding top x-axis) shows the resulting 1D radiative-convective equilibrium temperature profile.The dashed lines (and bottom x-axis) show the vertical gas mixing ratio profiles under the assumption of thermochemical equilibrium.These abundances, along with the absorption crosssections shown in the bottom panel of Figure 3, are what control the relative contributions of each gaseous opacity to the total transmission spectrum.

PICASO
The core 1D radiative-convective model is based upon the legacy "Extrasolar Giant Planet (EGP)" code described in Refs.[69, 80, 83] and since updated and modernized within the PICASO 84 framework described in Mukherjee et al., submitted (PICASO 3.0).The PICASO 3.0 model uses gaseous opacities created from the references listed in Ref [80].The grid of PICASO models contains metallicity points at 0.1, 0.3, 1, 3, 10, 30, 50, and 100x solar; C/O at 0.23, 0.46, 0.69, and 0.92; and also assumes full day-night heat redistribution.The clouds are modeled using the Virga 85 implementation of the Eddysed 86 framework, which requires a vertical mixing coefficient (constant with altitude; log10Kzz = 5, 7, 9, and 11 [cgs units]) and a vertically-constant sedimentation parameter (fsed = 0.6, 1, 3, 6, and 10), with optical/material properties for clouds thought to exist at WASP-39b's pressures and temperatures (Na2S, MnS, and MgSiO3).The fsed parameter controls the vertical extent of the cloud, and Kzz and fsed together control the mean droplet sizes with altitude in the atmosphere.A chi-square grid search along the described dimensions is performed to identify the best fit.Within this grid, the nominal best fit ( 2 / "#$# =1.34) is 10x solar metallicity, a sub-solar C/O (0.23), with an extended large droplet cloud (fsed = 0.6, log10Kzz = 9) that produces a grey continuum over these wavelengths, consistent with the ScCHIMERA results above.

ATMO
The ATMO radiative-convective-thermochemical equilibrium solver is described in Refs [87-90].This grid consists of model transmission spectra for four different day-night energy redistribution factors (0.25, 0.5, 0.75, 1.0, where 0.5 is "full", and 1.0 is "dayside-only"), six metallicities (0.1, 1, 10, 50, 100, 200 times solar), six C/O ratios (0.35, 0.55, 0.70, 0.75, 1.0, 1.5), two haze factors (no haze and 10 times multi-gas Rayleigh Scattering) and four grey cloud factors (no cloud, 0.5, 1, and 5 times the strength of H2 Rayleigh Scattering at 350 nm between 1 and 50 mbar pressure levels).Each model transmission spectrum from the grid is binned to the same resolution as that of the observations to compute  2 , with a (wavelength-independent) transit depth offset as the free parameter.Within this grid, we find a best-fit model ( 2 / "#$# = 1.39) spectrum arising from a redistribution factor of 0.75 (slightly hotter than a full day-night redistribution would produce), a metallicity of 10x solar, a super-solar C/O ratio of 0.7, a haze factor of 10, and a cloud factor of 5.

Quantifying Feature Detection Significance
We quantified the detection significance 95 of CO2 with the following steps: 1) the best-fit grid model without CO2 (i.e., the "no CO2" black curve shown in Figure 3) is first subtracted from the data, leaving behind a strong residual feature due to CO2 (ED Figure 3).The peak per-spectral-bin mean SNR of this residual feature is ~10σ.To utilize the full line/band shape we then fit the residual peak with (1) a four-parameter Gaussian model (centroid, amplitude, width, and vertical offset), shown as red curves in ED Figure 3, and (2) a "no feature" constant using a nested sampling routine 74 .The Bayesian evidence between the Gaussian model and constant model were then used to compute a Bayes factor, B, and corresponding detection significance 96 .For the CO2 residual feature, ln(B) is 340.5, which equates to a 26.2σ detection.From this analysis, we conclude that the CO2 feature is robustly detected.
Upon inspecting Figures 2 and 3 in the main text, there appears to be a feature near 4.0 µm (just short of the major CO2 feature).We repeated the same analysis as above, but instead compared the Bayesian evidence from a 2-component Gaussian model fit (to accommodate for both the CO2 feature and the unknown absorber) to that of the single component Gaussian model fit above.Upon doing so, we find ln(B) = 0.98 which equates to a 2σ significance.Restricting the prior range for the second Gaussian to be localized near the 4 µm feature boosts the significance to 2.3σ.Future analyses will focus on the nature of this feature and more rigorous quantification via nested Bayesian model comparison within atmospheric retrieval frameworks (e.g., Ref [34]).
ED Figure 3: Assessment of the strength of spectral features for WASP-39b.Residual features (blue data points) after subtracting the continuum best model (black "no CO2" model curve in Figure 3).A best-fitting ensemble of a 2-component Gaussian model to both the CO2 feature and the unknown absorber feature (~4 µm) is shown in red.
as noted in the text.All the data and models presented in this publication can be found at https://doi.10.5281/zenodo.6959427.the design of the program.KBS generated the observing plan with input from the team.EvSc, NE, and TGB provided instrument expertise.BB, EK, HW, IC, JLB, LK, MLM, MRL, NMB, and ZBT led or co-led working groups and/or contributed to significant strategic planning efforts like the design and implementation of the pre-launch Data Challenges.AC, DS, EvSc, NE, NG, TGP, VP generated simulated data for pre-launch testing of methods.AC, AF, CP, EA, EvSc, JaK, LA, TJB, and ZR contributed to the development of data analysis pipelines and/or provided the data analysis products used in this analysis.Reduced the data, modeled the light curves, and produced the planetary spectrum.JG, JoL, KO, MRL, NEB, SaMu, and SEM generated theoretical model grids for comparison with data.AC, DS, EvSc, HW, JaK, JF, JG, JLB, JoL, KBS, KO, MRL, NEB, NMB, SaMu, SEM, ZBT, and ZR contributed significantly to the writing of this manuscript.EA, EvSc, MRL, and ZBT generated figures for this manuscript.

Figure 1 :
Figure 1: JWST NIRSpec time-series data for WASP-39b.a) Spectroscopic light curves for WASP-39b's transit with a spectral resolving power of 20 and a time cadence of 1 minute (data are binned and offset vertically for display purposes only).An exoplanet light curve model was fitted to the data using a quadratic limb darkening law with an exponential ramp and a quadratic function of time removed.b) Residuals of the binned light curve after subtracting the transit model scaled up by a factor of five to show the structure.The RMS of the residuals are given in units of ppm.The numbers in brackets are the ratio of the RMS to the predicted photon-limited noise.

Figure 3 :
Figure 3: Interpretation of WASP-39b's transmission spectrum.The top panel shows a comparison of the FIREFLy reduction to the best-fit ScCHIMERA theoretical model binned to the resolution of the data (blue curve, see Methods).The key parameters of the model are 10x solar metallicity, carbon-to-oxygen ratio of 0.35, and cloud opacity of 7x10 -3 cm 2 /g.The impact of the opacity sources expected from thermochemical equilibrium over the full bandpass are indicated by removing the opacity contribution from individual gasses one at a time.As in Figure2, the axis on the right shows equivalent scale heights in WASP-39b's atmosphere.The bottom panel shows the molecular absorption cross-sections for each gas in the best-fit model.The model is well matched to the data ( !/ "#$# =1.3), suggesting that our assumptions broadly capture the important physics and chemistry in WASP-39b's atmosphere.However, there is a feature near 4.0 µm that cannot be reproduced by the models used here.The strong CO2 absorption (4.1 -4.6 µm) and the apparent lack of methane (3.0 -3.5 µm) is what drives the solution to an elevated atmospheric metal enrichment, ruling out previous low metallicity estimates[29][30][31] .The other reductions and models give similar results.