Measuring the abundances of carbon and oxygen in exoplanet atmospheres is considered a crucial avenue for unlocking the formation and evolution of exoplanetary systems^{1,2}. Access to the chemical inventory of an exoplanet requires high-precision observations, often inferred from individual molecular detections with low-resolution space-based^{3–5} and high-resolution ground-based^{6–8} facilities. Here we report the medium-resolution (^{9}), obtained with the Near Infrared Spectrograph (NIRSpec) G395H grating of JWST. Our observations achieve 1.46 times photon precision, providing an average transit depth uncertainty of 221 ppm per spectroscopic bin, and present minimal impacts from systematic effects. We detect significant absorption from CO_{2} (28.5_{2}O (21.5_{2} as the source of absorption at 4.1 μm (4.8_{2}, underscore the importance of characterizing the chemistry in exoplanet atmospheres and showcase NIRSpec G395H as an excellent mode for time-series observations over this critical wavelength range^{10}.

The medium-resolution transmission spectrum of the exoplanet WASP-39b, described using observations from the Near Infrared Spectrograph G395H grating aboard JWST, shows significant absorption from CO_{2} and H_{2}O and detection of SO_{2}.

We obtained a single-transit observation of WASP-39b using the NIRSpec^{11,12} G395H grating on 30–31 July 2022 between 21:45 and 06:21 UTC using the Bright Object Time Series mode. WASP-39b is a hot (_{eq} = 1,120 K), low-density giant planet with an extended atmosphere. Previous spectroscopic observations have shown prominent atmospheric absorption by Na, K and H_{2}O (refs. ^{3,4,13–15}), with tentative evidence of CO_{2} from infrared photometry^{4}. Atmospheric models fitted to the spectrum have inferred metallicities (amount of heavy elements relative to the host star) from 0.003 to 300 times solar^{3,15–20}, which makes it difficult to ascertain the formation pathway of the planet^{21,22}. The host, WASP-39, is a G8-type star that shows little photometric variability^{23} and has nearly solar elemental abundance patterns^{24}. The quiet host and extended planetary atmosphere make WASP-39b an ideal exoplanet for transmission spectroscopy^{25}. The transmission spectrum of WASP-39b was observed as part of the JWST Transiting Exoplanet Community Director’s Discretionary Early Release Science (JTEC ERS) Program^{26,27} (ERS-1366; principal investigators Natalie M. Batalha, Jacob L. Bean and Kevin B. Stevenson), which uses four instrument configurations to test their capabilities and provide lessons learned for the community.

The NIRSpec G395H data were recorded with the 1.6″ × 1.6″ fixed slit aperture using the SUB2048 subarray and NRSRAPID readout pattern, with spectra dispersed across both the NRS1 and NRS2 detectors. Over the roughly 8-h duration of the observation, a total of 465 integrations were taken, centred around the 2.8-h transit. We obtained 70 groups per integration, resulting in an effective integration time of 63.14 s. During the observation, the telescope experienced a ‘tilt event’, a spontaneous and abrupt change in the position of one or more mirror segments, causing changes in the point spread function (PSF) and hence jumps in flux^{28}. The tilt event occurred mid-transit, affecting approximately three integrations and resulted in a noticeable step in the flux time series, the size of which is dependent on wavelength (Fig.

We produced several reductions of the observations using independent analysis pipelines (see

We show transmission spectra from several combinations of independent reductions and light-curve-fitting routines in Fig.

We show the resultant spectra from five out of 11 independent fitting pipelines, which used distinct analysis methods to demonstrate the robust structure of the spectrum (see ^{2}, in which

_{2}O and CO_{2} with a grey-cloud-top pressure corresponding to ≈1 mbar. The models find that the data are best explained by 3–10 times solar metallicity (M/H) and sub-solar to solar C/O (C/O = 0.30–0.46). The extra absorption owing to SO_{2}, seen in the spectrum around 4.1 μm, is not included in the RCTE model grids and causes a marked impact on the ^{2}/

We compared the weighted-average G395H transmission spectrum to three grids of 1D radiative–convective–thermochemical equilibrium (RCTE) atmosphere models of WASP-39b (see ^{2}/_{2}O (at 21.5^{29}) can be attributed to CO_{2} (28.5_{2}O and CO_{2} are expected atmospheric constituents for near-solar atmospheric metallicities, with the CO_{2} abundance increasing nonlinearly with higher metallicity^{30}. The spectral feature at 4.56 μm (3.3^{31}, we interpret this feature as SO_{2} (4.8

Although SO_{2} would have volume mixing ratios (VMRs) of less than 10^{−10} throughout most of the observable atmosphere in thermochemical equilibrium, coupled photochemistry of H_{2}S and H_{2}O can produce SO_{2} on giant exoplanets, with the resulting SO_{2} mixing ratio expected to increase with increasing atmospheric metallicity^{32–34}. We find that a VMR of approximately 10^{−6} of SO_{2} is required to fit the spectral feature at 4.1 μm in the transmission spectrum of WASP-39b, consistent with lower-resolution NIRSpec PRISM observations of this planet^{31} and previous photochemical modelling of super-solar metallicity giant exoplanets^{34,35}. Figure ^{2}/^{−5.6} injected SO_{2}. The inclusion of SO_{2} in the models results in an improved ^{2}/

^{2}/^{−5.6} (VMR) SO_{2}. We also show this model with a selection of the anticipated absorbing species and the cloud opacity removed to indicate their contributions to the model. The inclusion of SO_{2} in the model decreases the ^{2}/_{2}S, OCS and CH_{4}, although their spectral features cannot be robustly detected in the spectrum. We show a model without CO and CH_{4} in

We also look for evidence of CH_{4}, CO, H_{2}S and OCS (carbonyl sulfide) because their near-solar chemical equilibrium abundances could result in a contribution to the spectrum. We see no evidence of CH_{4} in our spectrum between 3.0 and 3.6 μm (ref. ^{23}), which is indicative of C/O < 1 (ref. ^{36}) and/or photochemical destruction^{34,37}. With regards to CO, H_{2}S and OCS, we were unable to conclusively confirm their presence with these data. In particular, CO, H_{2}O, OCS and our modelled cloud deck all have overlapping opacity, which creates a pseudo-continuum from 4.6 to 5.1 μm (see Figs.

Our models show an atmosphere enriched in heavy elements, with best-fit parameters ranging from 3 to 10 times solar metallicity, given the spacing of individual model grids (see ^{38} but 10 times solar metallicity when assuming a grey cloud. In general, forward model grids fit the main features of the data but do not place statistically significant constraints on many of the atmospheric parameters (see

We are able to strongly rule out an absolute C/O ≥ 1 scenario (^{2}/_{2} ice line at which the gas may be carbon-rich^{39}. Our C/O upper limit, therefore, suggests that WASP-39b may have either formed at smaller orbital radii with gas-dominated accretion or that the accretion of solids enriched the atmosphere of WASP-39b with oxygen-bearing species^{2}. The level of metal enrichment (3–10 times solar) is consistent with similar measurements of Jupiter and Saturn^{40,41}, potentially suggesting core-accretion formation scenarios^{42}, and is consistent with upper limits from interior-structure modelling^{43}. These NIRSpec G395H transmission spectroscopy observations demonstrate the promise of robustly characterizing the atmospheric properties of exoplanets with JWST unburdened by substantial instrumental systematics, unravelling the nature and origins of exoplanetary systems.

We produced several analyses of stellar spectra from the Stage 1 2D spectral images produced using the default STScI JWST Calibration Pipeline^{44} (‘rateints’ files) and by means of customized runs of the STScI JWST Calibration Pipeline with user-defined inputs and processes for steps such as the ‘jump detection’ and ‘bias subtraction’ steps.

Each pipeline starts with the raw ‘uncal’ 2D images that contain group-level products. As we noticed that the default superbias images were of poor quality, we produced two customized runs of the JWST Calibration Pipeline, using either the default bias step or a customized version. The customized step built a pseudo-bias image by computing the median pixel value in the first group across all integrations and then subtracted the new bias image from all groups. We note that the poor quality of the default superbias images affects NRS1 more notably than NRS2, and this method could be revised once a better superbias is available.

Before ramp fitting, both our standard and custom bias step runs of the edited JWST Calibration Pipeline ‘destriped’ the group-level images to remove so-called ‘1/^{11,12}). Our group-level destriping step used a mask of the trace 15

The results of our customized runs of the JWST Calibration Pipeline are a set of custom group-level destriped products and custom bias-subtracted group-level destriped products. In both cases, the ramp-jump detection threshold of the JWST Calibration Pipeline was set to 15

For all analyses, wavelength maps from the JWST Calibration Pipeline were used to produce wavelength solutions, verified against stellar absorption lines, for both detectors. The mid-integration times in BJD_{TDB} were extracted from the image headers for use in producing light curves. None of our data-reduction pipelines performed a flat-field correction, as the available flat fields were of poor quality and unexpectedly removed portions of the spectral trace. In general, we found that 1/

Below we detail each of the independent data-reduction pipelines used to produce the time series of stellar spectra from our G395H observations.

We used the Exoplanet Timeseries Characterisation - JWST Extraction and Diagnostics Investigator (ExoTiC-JEDI^{45}) pipeline on our custom group-level destriped products, treating each detector separately. Using the data-quality flags produced by the JWST Calibration Pipeline, we replaced any pixels identified as bad, saturated, dead, hot, low quantum efficiency or no gain value with the median value of surrounding pixels. We also searched each integration for pixels that were spatial outliers from the median of the surrounding 20 pixels in the same row by 6^{−4} and 8 × 10^{−6} of a pixel for NRS1 and NRS2, respectively. The aperture column sums resulted in 1D stellar spectra with errors calculated from photon noise after converting from data numbers using the gain factor. This reduction is denoted hereafter as ExoTiC-JEDI [V1].

We produced further 1D stellar spectra from both the custom group-level destriped product and custom bias-subtracted group-level destriped products using the ExoTiC-JEDI pipeline as described above, but with further cleaning by repeating the spatial outliers step. The reduction with further cleaning using the custom group-level destriped products is hence denoted as ExoTiC-JEDI [V2] and the reduction with further cleaning using the custom bias-subtracted group-level destriped products is hence denoted as ExoTiC-JEDI [V3].

We used the Tiberius pipeline, which builds on the LRG-BEASTS spectral reduction and analysis pipelines^{15,46,47}, on our custom group-level destriped products. For each detector, we created bad-pixel masks by manually identifying hot pixels in the data. We then combined them with pixels flagged as greater than 3

We traced the spectra by fitting Gaussians at each column and used a running median, calculated with a moving box with a width of five data points, to smooth the measured centres of the trace. We fitted these smoothed centres with a fourth-order polynomial, removed points that deviated from the median by 3

We used the transitspectroscopy pipeline^{48} on the ‘rateints’ products of the JWST Calibration Pipeline, treating each detector separately. The trace position was found from the median integration by cross-correlating each column with a Gaussian function, removing outliers using a median filter with a 10-pixel-wide window and smoothing the trace with a spline. We removed 1/^{49} to obtain the 1D stellar spectra, with a 5-pixel-wide aperture above and below the trace. This allowed us to treat bad pixels and cosmic rays that had not been accounted for or masked in the ‘rateints’ products in an automated fashion. To monitor systematic trends in the observations, we also calculated the trace centre as described above and the FWHM for all integrations. The FWHM was calculated at each column and at each integration by first subtracting each column to half the maximum value in it, with a spline used to interpolate the profile. The roots of this profile were then found to estimate the FWHM.

We used two customized versions of the Eureka! pipeline^{50}, which combines standard steps from the JWST Calibration Pipeline with an optimal extraction scheme to obtain the time series of stellar spectra.

The first Eureka! reduction used our custom group-level destriped products and applied Stages 2 and 3 of Eureka! Stage 2, a wrapper of the JWST Calibration Pipeline, followed the default settings up to the flat fielding and photometric calibration steps, which were both skipped. Stage 3 of Eureka! was then used to perform the background subtraction and extraction of the 1D stellar spectra. We started by correcting for the curvature of NIRSpec G395H spectra by shifting the detector columns by whole pixels, to bring the peak of the distribution of the counts in each column to the centre of our subarray. To ensure that this curvature correction was smooth, we computed the shifts in each column for each integration from the median integration frame in each segment and applied a running median to the shifts obtained for each column. The pixel shifts were applied with periodic boundary conditions, such that pixels shifted upwards from the top of the subarray appeared at the bottom after the correction, ensuring no pixels were lost. We applied a column-by-column background subtraction by fitting and subtracting a flat line to each column of the curvature-corrected data frames, obtained by fitting all pixels further than six pixels from the central row. We also performed a double iteration of outlier rejection in time with a threshold of 10

The second Eureka! reduction (Eureka! [V2]) used the ‘rateints’ outputs of the JWST Calibration Pipeline and applied Stage 2 of Eureka! as described above, with a modified version of Stage 3. In this reduction, we corrected the curvature of the trace using a spline and found the centre of the trace using the median of each column. We removed 1/

Limb-darkening is a function of the physical structure of the star that results in variations in the specific intensity of the light from the centre of the star to the limb, such that the limb looks darker than the centre. This is because of the change in depth of the stellar atmosphere being observed. At the limb of the star, the region of the atmosphere being observed at slant geometry is at higher altitudes and lower density, and thus lower temperatures, compared with the deeper atmosphere observed at the centre of the star, at which hotter, denser layers are observed. The effect of limb-darkening is most prominent at shorter wavelengths, resulting in a more U-shaped light curve compared with the flat-bottomed light curves observed at longer wavelengths. To account for the effects of limb-darkening on the time-series light curves, we used analytical approximations for computing the ratio of the mean intensity to the central intensity of the star. The most commonly used limb-darkening laws for exoplanet transit light curves are the quadratic, square-root and nonlinear four-parameter laws^{51}:

Quadratic:

Square-root:

Nonlinear four-parameter:

in which _{1}, _{2}, _{1}, _{2}, _{1}, _{2}, _{3} and _{4} are the limb-darkening coefficients and

The limb-darkening coefficients depend on the particular stellar atmosphere and therefore vary from star to star. For consistency across all of the light-curve fitting, we used 3D stellar models^{52} for _{eff} = 5,512 K, log(^{53}, a custom throughput was produced from the median of the measured ExoTiC-JEDI [V2] stellar spectra, divided by the stellar model and Gaussian smoothed.

For the limb-darkening coefficients that were held fixed, we used the values computed using the ExoTiC-LD^{54,55} package, which can compute the linear, quadratic and three-parameter and four-parameter nonlinear limb-darkening coefficients^{51,56}. To compute and fit for the coefficients from the square-root law, we used previously outlined formalisms^{57,58}. We highlight that we do not see any dependence in our transmission spectra on the limb-darkening procedure used across our independent reductions and analyses.

From the time series of extracted 1D stellar spectra, we created our broadband transit light curves by summing the flux over 2.725–3.716 μm for NRS1 and 3.829–5.172 μm for NRS2. For the spectroscopic light curves, we used a common 10-pixel binning scheme within these wavelength ranges to generate a total of 349 spectroscopic bins (146 for NRS1 and 203 for NRS2). We also tested wider and narrower binning schemes but found that 10-pixel-wide bins achieved the best compromise between the noise in the spectrum and showcasing the abilities of G395H across analyses. In our analyses, we treated the NRS1 and NRS2 light curves independently to account for differing systematics across the two detectors. To construct the full NIRSpec G395H transmission spectrum of WASP-39b, we fitted the NRS1 and NRS2 broadband and spectroscopic light curves using 11 independent light-curve-fitting codes, which are detailed below. When starting values were required, all analyses used the same system parameters^{37}. In many of our analyses, we detrended the raw broadband and spectroscopic light curves using the time-dependent decorrelation parameters for the change in the FWHM of the spectral trace or the shift in

Using fitting pipeline 1, we measured a centre of transit time (_{0}) of _{0} = 2,459,791.612039 ± 0.000017 BJD_{TDB} and _{0} = 2,459,791.6120689 ± 0.000021 BJD_{TDB} computed from the NRS1 and NRS2 broadband light curves, respectively; most of the fitting pipelines obtained _{0} within 1

For each of our analyses, we computed the expected photon noise from the raw counts taking into account the instrument read noise (16.18 e^{−} on NRS1 and 17.75 e^{−} on NRS2), gain (1.42 for NRS1 and 1.62 for NRS2) and the background counts (which are found to be negligible after cleaning) and compared it to the final signal-to-noise ratio in our light curves (see Fig. ^{59}, which is used to measure the deviation from the expected photon noise by binning the data into successively smaller bins (that is, fewer data points per bin) and calculating the signal-to-noise ratio achieved^{60}. Extended Data Fig. ^{54}).

Although there is a general consensus across each of the data analyses, by comparing the results of each fitting pipeline, we were better able to evaluate the impact of different approaches to the data reduction, such as the removal of bad pixels. For future studies, we recommend the application of several pipelines that use differing analysis methods, such as the treatment of limb-darkening, systematic effects and noise removal. No single pipeline presented on its own can fully evaluate the measured impact of each effect, given the differing strategies, targets and potential for chance events such as a mirror tilt with each observation. In particular, attention should be paid to 1/

Below, we detail each of our 11 fitting pipelines and summarise them in Extended Data Table

We fitted the broadband and spectroscopic light curves produced from the ExoTIC-JEDI [V3] stellar spectra using the least-squares optimizer, scipy.optimize lm (ref. ^{61}). We simultaneously fitted a batman transit model^{62} with a constant baseline and systematics models for data pre-tilt and post-tilt event, fixing the centre of transit time _{0}, the ratio of the semi-major axis to stellar radius _{⋆} and the inclination

We used the broadband light curves generated from the Tiberius stellar spectra and fitted for the ratio of the planet to stellar radii _{p}/_{⋆}, as well as _{0}, _{⋆}, the quadratic law limb-darkening coefficient _{1} and the systematics model parameters, the _{2} fixed. We used uniform priors for all the fitted parameters. Our analytic transit light-curve model was generated with batman. We fitted our broadband light curve with a transit + systematics model using a Gaussian process (GP)^{63,64}, implemented through george^{65}, and a Markov chain Monte Carlo method, implemented through emcee^{66}. For our Tiberius spectroscopic light curves, we held _{⋆}, _{0} fixed to the values determined from the broadband light-curve fits and applied a systematics correction from the broadband light-curve fit to aid in fitting the mirror-tilt event. We fitted the spectroscopic light curves using a GP with an exponential squared kernel for the same systematics detrending parameters detailed above. We used a Gaussian prior for _{⋆} and uniform priors for all other fitted parameters.

We used transit light curves from the ExoTiC-JEDI [V1] stellar spectra and fit the broadband and spectroscopic light curves using the least-squares minimizer LMFIT^{67}. We fitted each light curve with a two-component function consisting of a transit model (generated using batman) multiplied by a systematics model. Our systematics model included the ^{2} and ^{2}). To capture the amplitude of the tilt event in our systematics model, we carried out piecewise linear regression on the out-of-transit baseline pre-tilt and post-tilt. We first fit the broadband light curve by fixing _{0}, _{⋆}, _{p}/_{⋆}, stellar baseline flux and systematic trends using wide uniform priors. For the spectroscopic light curves, we fixed _{0}, _{⋆} and _{p}/_{⋆}. We held the nonlinear limb-darkening coefficients fixed.

We fit the broadband and spectroscopic light curves produced from the transitspectroscopy stellar spectra, running juliet^{68} in parallel with the light-curve-fitting module of the transitspectroscopy pipeline^{48} with dynamic nested sampling through dynesty^{69} and analytical transit models computed using batman. We fit the broadband light curves for NRS1 and NRS2 individually, fixing _{0}, _{⋆}, _{p}/_{⋆}, extra jitter and the mean out-of-transit flux. We also fitted two linear regressors, a simple slope and a ‘jump’ (a regressor with zeros before the tilt event and ones after the tilt event), scaled to fit the data. We fitted the square-root-law limb-darkening coefficients using the Kipping sampling scheme. We fitted the spectroscopic light curves at the native resolution of the instrument, fixing _{0}, _{⋆} and ^{70} and passing them through the SPAM algorithm^{71}. We binned the data into 10-pixel-wavelength bins after fitting the native-resolution light curves.

We fitted the transit light curves from the Eureka! [V1] stellar spectra using the ExoTEP analysis framework^{72–75}. ExoTEP uses batman to generate analytical light-curve models, adds an analytical instrument systematics model along with a photometric scatter parameter and fits for the best-fit parameters and their uncertainties using emcee. Before fitting, we cleaned the light curves by running ten iterations of 5_{p}/_{⋆}, photometric scatter, _{0}, _{⋆}, the quadratic limb-darkening coefficients and the systematics model parameters (normalization constant, slope in time and constant ‘jump’ offset). We used uninformative flat priors on all the parameters. The orbital parameters were fixed to the best-fit broadband light curve values for the subsequent spectroscopic light-curve fits.

We fitted transit light curves from the ExoTiC-JEDI [V1] stellar spectra using a custom lmfit light-curve-fitting code. The final systematic detrending model included a batman analytical transit model multiplied by a systematics model consisting of a linear stellar baseline term, a linear term for the _{0}, _{p}/_{⋆}, _{⋆}, _{0}, _{⋆} and the exponential ramp timescale to the broadband light-curve-fit values, and fitted for _{p}/_{⋆}, the

We fitted transit light curves from the Eureka! [V2] stellar spectra, using PyLightcurve (ref. ^{75}) to generate the transit model with emcee as the sampler. We calculated the nonlinear four-parameter limb-darkening coefficients using ExoTHETyS (ref. ^{76}), which relies on PHOENIX 2012–2013 stellar models^{77,78}, and fixed these in our fits to the precomputed theoretical values. Our full transit + systematics model included a transit model multiplied by a second-order polynomial in the time domain. We accounted for the tilt event by subtracting the mean of the last 30 integrations of the pre-transit data from the mean of the first 30 integrations of the post-transit data, to account for the jump in flux, shifting the post-transit light curve upwards by the jump value. We fitted for the systematics (the parameters of the second-order polynomial), _{p}/_{⋆} and _{0}. We used uniform priors for all the fitted parameters. We adopted the root mean square of the out-of-transit data as the error bars for the light-curve data points to account for the scatter in the data.

We used the transit light curves generated from the ExoTiC-JEDI [V1] stellar spectra. We fit the broadband light curves with a batman transit model multiplied by a second-order systematics model as a function of _{p}/_{⋆}, _{0} and _{⋆}, using wide uninformed priors, and ran our fits using emcee. For the spectroscopic light-curve fits, we fixed _{⋆} to the broadband light-curve best-fit values and fitted for an extra error term added in quadrature.

We used the transit light curves from the ExoTiC-JEDI [V1] stellar spectra. We fixed both of the quadratic limb-darkening coefficients and fitted the light curves with a batman transit model multiplied by a systematics model of a second-order function of _{0}, _{⋆} and _{p}/_{⋆} using emcee for each 10-pixel bin, with the walkers initialized in a tight cluster around the best-fit solution from a Levenberg–Marquardt minimization. For both the broadband and spectroscopic light curves, we also fit for an extra per-point error term.

We fitted the transit light curves from the ExoTiC-JEDI [V2] stellar spectra and performed our model fitting using automatic differentiation implemented with JAX (ref. ^{79}). We used a GP systematics model with a time-dependent Matérn (^{80}, making use of previously developed light-curve models^{81,82}. For the GP systematics component, a generalization of the algorithm used by the celerite package^{83} was adapted for JAX. We fixed both of the quadratic limb-darkening coefficients. For the initial broadband light-curve fit, both NRS1 and NRS2 were fitted simultaneously. All transit parameters were shared across both light curves, except for _{p}/_{⋆}, which was allowed to vary for NRS1 and NRS2 independently. We fitted for _{0}, the transit duration _{p}/_{⋆} values. For the spectroscopic light-curve fits, all transit parameters were then fixed to the maximum-likelihood values determined from the broadband fit, except for _{p}/_{⋆}, which was allowed to vary for each wavelength bin. Uncertainties for the transit model parameters, including _{p}/_{⋆}, were assumed to be Gaussian and estimated using the Fisher information matrix at the location of the maximum-likelihood solution, which was computed exactly using the JAX automatic differentiation framework.

We used transit light curves from the Eureka! [V2] time-series stellar spectra with the open-source Eureka! code to estimate the best-fit transit parameters and their uncertainties using a Markov chain Monte Carlo method fit to the data (implemented by emcee). A linear trend in time was used as a systematics model and a step function was used to account for the tilt event. We fixed _{⋆}, _{0} and the time of the tilt event to the best-fit values from the NRS1 broadband light curve, with the three integrations during the tilt event clipped from the light curve. We fitted for _{p}/_{⋆}, both quadratic limb-darkening coefficients, the linear time trend and the magnitude of the step from the tilt event, with uniform priors for both the magnitude of the step and the limb-darkening coefficients.

On the basis of the independent light-curve fits described above, we produced 11 transmission spectra from our NIRSpec G395H observations using several analyses and fitting methods. Extended Data Table

We computed the standard deviation of the 11 spectra in each wavelength bin and compared this to the mean uncertainty obtained in that bin. The average standard deviation in each bin across all fitting pipelines was 199 ppm, compared with an average uncertainty of 221 ppm (which ranged from 131 to 625 ppm across the bins). The computed standard deviation in each bin across all pipelines ranged from 85 to 1,040 ppm, with greater than 98% of the bins having a standard deviation lower than 500 ppm. We see an increase in scatter at longer wavelengths, with the structure of the scatter following closely with the measured stellar flux, for which throughput beyond 3.8 μm combines with decreasing stellar flux. The unweighted mean uncertainty of all 11 transmission spectra follows a similar structure to the standard deviation, with the uncertainty increasing at longer wavelengths. The uncertainties from each fitting pipeline are consistent to within 3

From all 11 transmission spectra, we computed a weighted-average transmission spectrum using the transit depth values from all reductions in each bin weighted by 1/variance (1/^{2}, in which

We calculated normalized transmission spectrum residuals for each fitting pipeline by subtracting the weighted-average spectrum and dividing by the uncertainty in each bin. We generated histograms of the normalized transmission spectrum residuals and used the mean and standard deviation of the residuals to compute a normalized probability density function (PDF). We performed a Kolmogorov–Smirnov test on each of the normalized residuals and found that they are all approximately symmetric around their means, with normal distributions. This confirms that they are Gaussian in shape, with the null hypothesis that they are not Gaussian strongly rejected in the majority of cases (see Extended Data Fig.

The PDFs of the residuals indicate three distinct clusters of computed spectra based on their deviations from the mean and their spreads. The first cluster is negatively offset by less than 200 ppm and corresponds to fitting pipelines that used extracted stellar spectra and that underwent further cleaning steps (for example, ExoTiC-JEDI [V3]). The second cluster is positively offset from the mean by about 120 ppm and contains most of the transmission spectra produced. We see no obvious trends in this group to any specific reduction or fitting process. The final cluster is centred around the mean but has a broad distribution, suggesting a larger scatter both above and below the average transmission spectrum. This is probably the result of uncleaned outliers or marginal offsets between NRS1 and NRS2. These transmission spectra demonstrate that the 11 independent fitting pipelines are able to accurately reproduce the same transmission spectral feature structures, further highlighting the minimized impact of systematics on the time-series light curves. We suspect that the minor differences resulting from different reduction products and fitting pipelines are linked to the superbias and treatment of 1/

To identify spectral absorption features, we compared the resulting weighted-average transmission spectrum of WASP-39b to several 1D RCTE atmosphere models from three independent model grids. Each forward model is computed on a set of common physical parameters (for example, metallicity, C/O ratio, internal temperature and heat redistribution), shown in Extended Data Table ^{84}. Each model transmission spectrum from the grids was binned to the same resolution as that of the observations to compute the ^{2} per data point, with a wavelength-independent transit depth offset as the free parameter. In general, the forward model grids fit the main features of the data but are unable to place statistically significant constraints on many of the atmospheric parameters, owing to both the finite nature of the forward model grid spacing^{13} and the insensitivity of some of these parameters to the 3–5-μm transmission spectrum of WASP-39b (for example, >100 K differences in interior temperature provided nearly identical ^{2}/

We used the ATMO RCTE grid^{85–88}, which consists of model transmission spectra for different day–night energy redistribution factors, atmospheric metallicities, C/O ratios, haze factors and grey cloud factors with a range of line lists and pressure-broadening sources^{88}. In total, there were 5,160 models. Within this grid, we find the best-fit model to have 3 times solar metallicity, with a C/O ratio of 0.35 and a grey cloud opacity 5 times the strength of H_{2} Rayleigh scattering at 350 nm and a ^{2}/

We calculated a grid of transmission spectra using the PHOENIX atmosphere model^{89–91}, varying the heat redistribution of the planet, atmospheric metallicity, C/O ratio, internal temperature, the presence of aerosols and the atmospheric chemistry (equilibrium or rainout). Opacities used include the BT2 H_{2}O line list^{92}, as well as HITRAN for 129 other main molecular absorbers^{93} and Kurucz and Bell data for atomic species^{94}. The HITRAN line lists available in this version of PHOENIX are often complete only at room temperature, which may be the cause of the apparent shift in the CO_{2} spectral feature compared with the other grids that primarily use HITEMP and ExoMol lists. This shift is the cause of the difference in ^{2} between PHOENIX and the other model grids. In total, there were 1,116 models. Within this grid, the best-fit model has 10 times solar metallicity, a C/O ratio of 0.3, an internal temperature of 400 K, rainout chemistry and a cloud deck top at 0.3 mbar. The best-fit model has a ^{2}/

We used the open-source radiative–convective equilibrium model PICASO 3.0 (refs. ^{95,96}), which has its heritage in the Fortran-based EGP mode^{97,98}, to compute a grid of 1D pressure–temperature models for WASP-39b. The opacity sources included in PICASO 3.0 are listed in Extended Data Table _{2}O (ref. ^{99}), CO_{2} (ref. ^{100}), CH_{4} (ref. ^{101}) and CO (ref. ^{102}). The parameters varied in this grid of models include the interior temperature of the planet (_{int}), atmospheric metallicity, C/O ratio and the dayside-to-nightside heat redistribution factor (see Extended Data Table ^{98,103}. In total, there were 192 cloud-free models. We include the effect of clouds in two ways. First, we post-processed the pressure–temperature profile using the cloud model Virga^{95,104}, which follows from previously developed methodologies^{38}, in which we included three condensable species (MnS, Na_{2}S and MgSiO_{3}). Virga requires a vertical mixing parameter, _{zz} (cm^{2} s^{−1}), and a vertically constant sedimentation efficiency parameter, _{sed}. In general, _{sed} controls the vertical extent of the cloud opacity, with low values (_{sed} < 1) creating large, vertically extended cloud decks with small particle sizes. In total, there were 3,840 cloudy models. The best fit from our grid with Virga-computed clouds has 3 times solar metallicity, solar C/O (0.458) and _{sed} = 0.6, which results in a ^{2}/

As well as the grid fit, we also use the PICASO framework to quantify the feature-detection significance. In this method, we are able to incorporate clouds on the fly using the fitting routine PyMultiNest^{105}. We fit for each of the grid parameters using a nearest-neighbour technique and a radius scaling to account for the unknown reference pressure, giving five parameters in total. When fitting for clouds, we either fit for _{zz} and _{sed} in the Virga framework (seven parameters in total) or we fit for the cloud-top pressure corresponding to a grey cloud deck with infinite opacity (six parameters in total). These results are described in the following section.

From the chemical equilibrium results of the single best-fit models, the molecules that could potentially contribute to the spectrum based on their abundances and 3–5-μm opacity sources are H_{2} and He (via continuum) and CO, H_{2}O, H_{2}S, CO_{2} and CH_{4}. More minor sources of opacity with VMR abundances <1 ppm are molecules such as OCS and NH_{3}. For example, removing H_{2}S, NH_{3} and OCS from the single best-fit PICASO 3.0 model increases the chi-square value by less than 0.002. Therefore, we focus on computing the statistical significance of only H_{2}O, SO_{2}, CO_{2}, CH_{4} and CO.

To quantify the statistical significance, we performed two different tests. First, we used a Gaussian residual fitting analysis, as used in other JTEC ERS analyses^{23,29,31}. In this method, we subtracted the best-fit model without a specific opacity source from the weighted-average spectrum of WASP-39b, isolating the supposed spectral feature. We then fit a three-parameter or four-parameter Gaussian curve to the residual data using a nested sampling algorithm to calculate the Bayesian evidence^{106}. For H_{2}O and CO, the extra transit depth offset parameter for the Gaussian fit was necessary to account for local mismatch of the fit to the continuum, whereas only a mean, standard deviation and scale parameter were required for a residual fit to the other molecules. We then compared this to the Bayesian evidence of a flat line to find the Bayes factor between a model that fits the spectral feature versus a model that excludes the spectral feature. These fits are shown in Extended Data Fig.

Although the Gaussian residual fitting method is useful for quantifying the presence of potentially unknown spectral features, it cannot robustly determine the source of any given opacity. We therefore used the Bayesian fitting routine from PyMultiNest in the PICASO 3.0 framework to refit the grid parameters, while excluding the opacity contribution from the species in question. Then, we compared the significance of the molecule through a Bayes factor analysis^{107}. Those values are shown in Extended Data Table

We find significant evidence (>3_{2}O, CO_{2} and SO_{2}. In general, the two methods only agree well for molecules whose contribution has a Gaussian shape (that is, SO_{2} and CO_{2}). For example, for CO_{2}, we find decisive 28.5_{2}O, we find 21.5_{2} is less substantial, but both methods give significant detections of 4.8^{31}. CO remained undetected with the Bayesian fitting analysis and therefore we are unable to robustly confirm evidence of CO. Similarly, no evidence for CH_{4} was found^{23}. Gaussian residual fitting in the region of CH_{4} absorption only found a very broad inverse Gaussian and so is not included in Extended Data Table

We performed an injection test with the PICASO best-fit model in the PyMultiNest fitting framework to determine the abundance of SO_{2} required to match the observations. We add SO_{2} opacity using the ExoMol line list^{108}, without rerunning the RCTE model to self-consistently compute a new climate profile. Fitting for the cloud deck dynamically, without SO_{2}, produces a single best estimate of 10 times solar metallicity, sub-solar C/O (0.229), resulting in a marginally worse ^{2}/_{2}, the single best fit tends back to 3 times solar metallicity, solar C/O. This suggests that cloud treatment and the exclusion of spectrally active molecules have an effect on the resultant physical interpretation of bulk atmospheric parameters. Ultimately, if we fit for SO_{2} in our PyMultiNest framework with the Virga cloud treatment, we obtain 3 times solar metallicity, solar C/O, log SO_{2} = −5.6 ± 0.1 (SO_{2} = 2.5 ± 0.65 ppm) and ^{2}/_{2}S, the presumed source for photochemically produced SO_{2} (ref. ^{36}).

To confirm the plausibility of SO_{2} absorption to explain the 4.1-μm spectral feature, we also computed models with prescribed, vertically uniform SO_{2} VMRs of 0, 1, 5 and 10 ppm using the structure from the best-fit PHOENIX model (10 times solar metallicity, C/O = 0.3). We calculated ad hoc spectra using the gCMCRT radiative transfer code^{109} with the ExoMol SO_{2} line list^{108} (see Extended Data Fig. _{2} abundance and performing a Levenberg–Marquardt regression gave a best-fit value of 4.6 ± 0.67 ppm. Inserting this abundance of SO_{2} into the best-fit PHOENIX model improves the ^{2}/

Future atmospheric retrievals can provide a more statistically robust measurement for the SO_{2} abundance and add extra information from the similar absorption seen in the PRISM transmission spectrum^{29,31}.

A 0.08-μm-wide bump in transit depth centred at 4.56 μm is not fit by any of the model grids. This feature, picked up by the resolution of G395H, is not clearly seen in other ERS observations of WASP-39b. Following the same Gaussian residual fitting procedure as described above, we found a feature significance of 3.3_{4} (ref. ^{110}), C_{2}H_{2} (ref. ^{111}), C_{2}H_{4} (ref. ^{112}), C_{2}H_{6} (ref. ^{113}), CO (ref. ^{114}), CO_{2} (ref. ^{100}), CS_{2} (ref. ^{113}), CN (ref. ^{115}), HCN (ref. ^{116}), HCl (ref. ^{113}), H_{2}S (ref. ^{117}), HF (ref. ^{118}), H_{3}^{+} (ref. ^{119}), LiCl (ref. ^{115}), NH_{3} (ref. ^{120}), NO (ref. ^{114}), NO_{2} (ref. ^{113}), N_{2}O (ref. ^{114}), N_{2} (ref. ^{121}), NaCl (ref. ^{122}), OCS (ref. ^{113}), PH_{3} (ref. ^{123}), PN (ref. ^{124}), PO (ref. ^{125}), SH (ref. ^{126}), SiS (ref. ^{127}), SiH_{4} (ref. ^{128}), SiO (ref. ^{129}), the X–X state of SO (ref. ^{130}), SO_{2} (ref. ^{108}), SO_{3} (ref. ^{108}) and isotopologues of H_{2}O, CH_{4}, CO_{2} and CO, but did not find a convincing candidate that showed opacity at the correct wavelength or the correct width. The narrowness of the feature suggests that it could be a very distinct Q-branch, in which the rotational quantum number in the ground state is the same as the rotational quantum number in the excited state. However, of the molecules we explored, there were no candidates with a distinct Q-branch at this wavelength whose P-branch and R-branch did not obstruct the neighbouring CO_{2} and continuum-like CO + H_{2}O opacity.

We also note that many of these species lack high-temperature line-list data, making it difficult to definitively rule out such species. For example, OCS, SO and CS_{2} are available in HITRAN2020 (ref. ^{113}) but not in ExoMol^{131}. Furthermore, if photochemistry is important for WASP-39b, as indicated by the presence of SO_{2}, then there may be many species out of equilibrium that may contribute to the transit spectrum, some of which do not have high-temperature opacity data at present (such as OCS, NH_{2} or HSO). Future observations over this wavelength region of this and other planets may confirm or refute the presence of this unknown absorber.

Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at

This work is based on observations made with the NASA/ESA/CSA JWST. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for the JWST. These observations are associated with programme JWST-ERS-01366. Support for programme JWST-ERS-01366 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127. L.A. acknowledges funding from STFC grant ST/W507337/1 and from the University of Bristol School of Physics PhD Scholarship Fund.

All authors played a substantial role in one or more of the following: development of the original proposal, management of the project, definition of the target list and observation plan, analysis of the data, theoretical modelling and preparation of this manuscript. Some specific contributions are listed as follows. N.M.B., J.L.B. and K.B.S. provided overall programme leadership and management. L.A. and H.R.W. led the efforts for this manuscript. D.K.S., E.M.-R.K., H.R.W., I.J.M.C., J.L.B., K.B.S., L.K., M.L.-M., M.R.L., N.M.B., V.P. and Z.K.B.-T. made notable contributions to the design of the programme. K.B.S. generated the observing plan, with input from the team. E.S. and N.E. provided instrument expertise. B.B., E.M.-R.K., H.R.W., I.J.M.C., J.L.B., L.K., M.L.-M., M.R.L., N.M.B. and Z.K.B.-T. led or co-led working groups and/or contributed to important strategic planning efforts, such as the design and implementation of the prelaunch Data Challenges. A.L.C., D.K.S., E.S., N.E., N.P.G. and V.P. generated simulated data for prelaunch testing of methods. L.A., H.R.W., M.K.A., N.E.B. and J.D.L. contributed substantially to the writing of this manuscript, along with contributions in

The data used in this paper are associated with JWST programme ERS 1366 (observation #4) and are available from the Mikulski Archive for Space Telescopes (MAST;

The codes used in this publication to extract, reduce and analyse the data are as follows; STScI JWST Calibration Pipeline^{44} (^{50} (^{45} (^{68} (^{15,46,47}, transitspectroscopy^{48} (^{62} (^{83} (^{69} (^{66} (^{80} (^{72–74}, ExoTHETyS^{76} (^{54} (^{55} (^{65} (^{79} (^{67} (^{75} (^{132}) (^{81} (^{133,134}, matplotlib^{135}, numpy^{136}, pandas^{137}, scipy^{61} and xarray^{138}. The atmospheric models used to fit the data can be found at ATMO^{85–88}, PHOENIX^{89–91}, PICASO^{95,96} (^{95,104} (^{109} (

The authors declare no competing interests.

_{2}O depletion found on

_{2}O abundances and cloud properties in ten hot giant exoplanets

_{2}O, Na, and K

_{eff}≤ 50000 K at several surface gravities

_{eff}≤ 4800 K

_{eff}≤ 10 000 K

_{eff}≤ 10,000 K

^{−}opacity as a probe in ultrahot Jupiters

_{2}isotopologues up to

^{−1}and 1500 K, with line shape parameters

_{4}and hot methane continuum

^{1}Σ

^{+}ground electronic state

_{2}

_{2}H

_{4}

_{2}S

_{3}

^{+}

_{2}from 4,500 to 15,700 cm

^{−1}revisited with PGOPHER

^{2}Π –

^{2}Π and

^{2}Σ

^{+}–

^{2}Π transitions of SH

^{139}.

^{139}.

Purple denotes NRS1 and orange denotes NRS2. The normalized flux offset is calculated per pixel by measuring the median flux in the stellar baseline before and after the transit and calculating the difference. These differences are then normalized by the before-transit flux and plotted on a common scale. Overplotted are the data binned to a resolution of 10 pixels to match our presented transmission spectra (Fig.

In each subplot, the red line shows the expected relationship for perfect Gaussian white noise. The black lines show the observed noise from each spectroscopic light curve for pipelines 1, 3 and 5. To compare bins and noise levels, values for all bins in each pipeline are normalized by dividing by the value for a bin width of 1.

Shown after all other absorption from the best-fit model is subtracted from the data. Each of the Gaussian fits has a higher Bayesian evidence than the flat-line fits, indicating a detection, although to varying degrees of significance.

Model transmission spectra compared with the observed spectral feature at 4.1 μm in the G395H data. At wavelengths short of 3.95 μm, which is outside the SO_{2} band, all models overlap, further suggesting that the data can be explained by the presence of SO_{2} in the atmosphere. By interpolating these 10 times solar metallicity models, we find a best-fit SO_{2} abundance of 4.6 ± 0.67 ppm. With the best-fit PICASO 3.0 at 3 times solar metallicity, we find an SO_{2} abundance of 2.5 ± 0.65 ppm.

Summary of transit light-curve fitting

Summary of transit light-curve fitting

An outline of the combined products and fitting pipelines used to compute each transmission spectrum.

RCTE model grids

RCTE model grids

The parameter space explored by each RCTE model grid. The best-fit model for each grid is shown in bold. Ref. ^{140–142}.

Detection significances

Detection significances

Feature-detection significance for dominant sources of opacity with two different methods.

Peer Review File

is available for this paper at

The online version contains supplementary material available at