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

The Saturn-mass exoplanet WASP-39b has been the subject of extensive efforts to determine its atmospheric properties using transmission spectroscopy1–4. However, these efforts have been hampered by modelling degeneracies between composition and cloud properties that are caused by limited data quality5–9. Here we present the transmission spectrum of WASP-39b obtained using the Single-Object Slitless Spectroscopy (SOSS) mode of the Near Infrared Imager and Slitless Spectrograph (NIRISS) instrument on the JWST. This spectrum spans 0.6–2.8 μm in wavelength and shows several water-absorption bands, the potassium resonance doublet and signatures of clouds. The precision and broad wavelength coverage of NIRISS/SOSS allows us to break model degeneracies between cloud properties and the atmospheric composition of WASP-39b, favouring a heavy-element enhancement (‘metallicity’) of about 10–30 times the solar value, a sub-solar carbon-to-oxygen (C/O) ratio and a solar-to-super-solar potassium-to-oxygen (K/O) ratio. The observations are also best explained by wavelength-dependent, non-grey clouds with inhomogeneous coverageof the planet’s terminator.

Transmission spectroscopy provides insight into the atmospheric properties and consequently the formation history, physics, and chemistry of transiting exoplanets 1 .However, obtaining precise inferences of atmospheric properties from transmission spectra requires simultaneously measuring the strength and shape of multiple spectral absorption features from a wide range of chemical species [2][3][4] .This has been challenging given the precision and wavelength coverage of previous observatories 5 .Here, we present the transmission spectrum of the Saturn-mass exoplanet WASP-39 b obtained using the SOSS mode of the NIRISS instrument on the JWST.This spectrum spans 0.6-2.8m in µ wavelength and reveals multiple water absorption bands, the potassium resonance doublet, as well as signatures of clouds.The precision and broad wavelength coverage of NIRISS-SOSS allows us to break model degeneracies between cloud properties and the atmospheric composition of WASP-39 b, favouring a heavy element enhancement ("metallicity") of ~10 -30 the solar value, a sub-solar carbon-to-oxygen (C/O) ratio, and × a solar-to-super-solar potassium-to-oxygen (K/O) ratio.The observations are also best explained by wavelength-dependent, non-gray clouds with inhomogeneous coverage of the planet's terminator.
We observed a transit of WASP-39 b using the Near Infrared Imager and Slitless Spectrograph (NIRISS) 6 on the JWST as part of the Transiting Exoplanet Community Early Release Science Program (ERS) 7,8 .Our observations spanned 8.2 hours starting on UTC July 26, 2022 20:45, covering the 2.8-hour transit as well as 3.0 hours prior to and 2.4 hours after the transit to establish a flux baseline.The data were taken in the Single Object Slitless Spectroscopy (SOSS) mode, which simultaneously covers the wavelength range from m across two 0. 6 − 2. 8 µ spectral orders on the same detector.Order 1 contains the spectral range between 0. 85 − 2. 8 µ m at an average resolving power of , while Order 2 delivers the spectral range  ≡ λ/∆λ = 700 of m at an average resolving power of .In the SOSS mode, the spectra 0. 6 − 1. 4 µ  = 1400 are spread across more than 20 pixels in the cross-dispersion direction via a cylindrical defocusing lens (see Extended Data Fig. 1), thus allowing longer integration times and reducing the impact of pixel-level differences in the detector response.However, this defocus also results in the physical overlap of both orders on the detector, as well as the characteristic "horned" structure of the SOSS PSF.The time series observation was composed of 537 integrations of 49.4 seconds (nine groups per integration), corresponding to a duty cycle of 89%.
We extracted the stellar spectra from the time series observations using six different pipelines to test the impact of differences in spectral order tracing, 1/ noise correction, background removal,  and spectrum extraction methodology (see Methods and Extended Data Figs. 2 and 3).For each pipeline, we created spectrophotometric light curves at either pixel resolution (that is one light curve per detector column) or the native instrument resolution (two to three pixels per bin) (Fig. 1).We also summed the data to create individual white-light curves for each order (Extended Data Fig. 4).The spectrophotometric and white-light curves are largely free of instrumental systematics except for a constant-rate linear trend in time and an exponential ramp effect within the first 15 minutes of the time series, the latter of which we trimmed from the remainder of the data analysis.The signal-to-noise in the spectrophotometric light curves for Order 1 at m is 165; the signal-to-noise for Order 2 at m is 103. 1. 34µ 0. 71 µ An exoplanet transit model (solid line) was fitted to each light curve with chromatic_fitting using a quadratic limb-darkening law.The limb-darkening coefficients, planet-to-star radius ratio (R p /R * ), and out-of-transit flux were varied in each wavelength channel, while all other parameters were fixed.The residuals to the best-fit models are shown for each light curve.The wavelength range per each channel is denoted on the left plot, while parts-per-million (ppm) scatter in the residuals is denoted on the right plot.We calculate the ppm as the standard deviation of the out-of-transit residuals.The reductions are from the nirHiss and chromatic_fitting routines described in the Methods. (https://github.com/afeinstein20/wasp39b_niriss_paper/blob/main/scripts/figure1.py) For each reduction (described in Methods), we fit the corresponding white-light curve per each order with a transit model to extract the best-fit orbital parameters (Extended Data Table 2).We fixed the orbital period to the best-fitting value from 9 .For the transit fits to the  spectrophotometric light curves, most fitting routines also set the time of mid-transit , impact  0 parameter , and semi-major axis to stellar radius ratio to the white-light parameters  / * (Extended Data Table 2) to ensure consistency between reductions.We thus allowed only the planet-to-star radius ratio, , and the two parameters of the quadratic limb-darkening law,   / *  1 and , to vary for each wavelength, along with any additional systematics and long-term trends  2 in the out-of-transit flux.These fits were carried out at either pixel resolution or the native instrument resolution depending on the pipeline.The high-resolution transit depths were subsequently binned into 100 spectral wavelength channels, 80 in Order 1 and 20 in Order 2, to create a transmission spectrum at .In Figure 2 we present this transmission spectrum of  ∼ 300 WASP-39b between m from the nirHiss, supreme-SPOON, and 0. 6 − 2. 8 µ transitspectroscopy reduction pipelines.We find consistent results between the different data reduction pipelines and light curve fitting analyses, with the derived spectra also in agreement with previous results from the Hubble Space Telescope (HST, see also Extended Data Fig. 5).
Figure 2: NIRISS transmission spectra for WASP-39b obtained by three data reduction pipelines.We find broad agreement in the overall structure of the transmission spectra between several reduction pipelines, a sample of which are presented here (see Extended Data Fig. 5 for all reductions).The JWST data are shown in the coloured points, while previous HST observations of WASP-39b 18 are shown in white.We note that we only consider wavelengths <0.85μm for Order 2, since Order 1 has much higher fidelity in the overlapping 0.85-1.0μmrange.The JWST and HST data agree across the three broad H 2 O features that they have in common.We find evidence of a K absorption feature at 0.76μm, which the new JWST data, which was hypothesised in the previous HST data 18 .(https://github.com/afeinstein20/wasp39b_niriss_paper/blob/main/scripts/figure2.py) We investigated the atmospheric properties of WASP-39b by comparing our measured transmission spectrum from the nirHiss pipeline to grids of one-dimensional, radiative-convective-thermochemical equilibrium models.These models explore the impact of atmospheric metallicity (M/H), carbon-to-oxygen ratio (C/O), potassium-to-oxygen ratio (K/O), heat redistribution (f), and cloud coverage on the transmission spectrum of the planet.We explored multiple cloud models ranging from parametric treatments 10,11 to a droplet sedimentation model 12 that calculates the vertical distributions of cloud mass mixing ratio and mean particle size from the balance between gravitational sedimentation and eddy diffusion of cloud particles.Using a Bayesian inference framework (see Methods), we compared these grids of models to the observations and inferred the range of M/H, C/O, K/O and f, that best explains the data while marginalising over different cloud treatments.WASP-39, the host star, has a metallicity equal to that of the Sun within measurement precision [13][14][15][16] , so we reference the planet's atmospheric abundances to the solar pattern of elemental abundances 17 .We compared the grid spectra computed by various models (PICASO, ATMO, Phoenix, and ScCHIMERA) with an observational spectrum obtained from each data reduction pipeline and obtained broadly consistent results on the inferred atmospheric properties.In what follows, we introduce the results from the comparison between the spectrum from nirHiss data reduction and ScCHIMERA grid that allows the most comprehensive treatments of cloud properties.
Our best-fitting model to the NIRISS-SOSS transmission spectrum of WASP-39 b is presented in Figure 3.This model comes from the ScCHIMERA grid (see Methods) and assumes an atmospheric metallicity of [M/H]=1.38 ( the solar value), C/O=0.2 ( solar value, 17 ), 23 × 0. 55 × [K/O]=0.1 (1.26 the solar value), full day-night heat redistribution (f=1), and flux-balanced × clouds with inhomogeneous terminator coverage.We label the unambiguous signatures of absorption due to potassium and water vapour in Figure 2. The spectral maxima at 0.9, 1.15, 1.4, and 1.8 m due to water absorption result in a detection of the molecule (see Methods).µ > 30σ Similarly, the potassium doublet at 0.768 m is detected in the data at .Besides these µ 6. 8σ chemical species, signatures of CO and/or CO 2 are identified due to their contribution to the spectrum past m, with a significance model preference for CO over CO 2 absorption 2. 3 µ 3. 6σ (see Methods).The two bottom panels show the molecular absorption cross-sections for a selection of gases observable within the NIRISS bandpass.The middle panel highlights gases inferred by our analysis of WASP-39b's spectrum.The bottom panel highlights some gases that were not identified in these data, but that may be present in future observations of other exoplanets.(https://github.com/afeinstein20/wasp39b_niriss_paper/blob/main/scripts/figure3.py) From the chemical equilibrium models considered, we find that the observations are best explained by a sub-solar C/O (see Fig. 4, panel a).Across the different spectroscopic resolutions and atmospheric models, the best-fit C/O is 0.2, which is the lowest ratio explored in the grid of models.We rule out super-solar C/O due to the lack of CH

5 − 5µ
We also find that the observations are best explained by an atmospheric metallicity of -10 30 × solar.Metallicity inferences over the wavelength range of these observations are largely driven by the size and shape of the water vapour features, with some minor contributions due to CO and/or CO 2 at longer wavelengths ( m; see Fig. 3 and Fig. 4, panel b).The preferred range > 2 µ of metallicities provides the best fit to the shape and size of the muted water vapour features shortward of m in combination with the larger water and CO/CO 2 feature longward of m, 2 µ 2 µ regardless of the assumed cloud treatment in our models.
Due to the simultaneous detection of potassium and water vapour, we are able to place constraints on the K/O ratio, which is a refractory-to-volatile elemental ratio, being a solar-to-super-solar value.Since the refractory elements are condensed into solids in most parts of protoplanetary disks, the disk gas accretion tends to cause a sub-stellar refractory elemental abundance 19 .By contrast, solid accretion, such as planetesimal accretion, acts to increase the refractory elemental abundance and refractory-to-volatile elemental ratio 20 , though the latter depends on the composition of the accreted solids 21 .We anticipate that the K/O ratio diagnoses to what degree the solid accretion enriched the atmosphere during the formation stage.All of our fitted models find that the WASP-39b observations are well described by solar-to-super-solar K/O ratios, which is in agreement with previous inferences for this planet obtained via observations with limited spectral coverage 22 .We do not expect the K feature to be impacted by stellar chromospheric magnetic activity given the effective temperature of the star K 23 ∼ 5300 and the general quietness of WASP-39 (see Ahrer et al. submitted) It is also in line with larger population studies of hot giant planets that broadly found solar-to-super-solar refractory abundances and solar-to-sub-solar H O abundances 22,24 .The shape and strength of the potassium  Under these equilibrium conditions, increasing the C/O results in less H 2 O and more CH 4 , the latter having spectroscopic signatures incompatible with the observations.To mute these incompatible CH 4 features at high C/O, the model requires a higher degree of cloudiness that also mutes any remaining H 2 O features in the spectrum.Bottom: The same exercise as above, but instead we vary the metallicity parameter.The metallicity constraint is driven by the λ > 2μm data; the high-metallicity models ([M/H]> 2) expect larger transit depths than what is seen in the data.The same reference model is plotted as a thick black line in both panels.(https://github.com/afeinstein20/wasp39b_niriss_paper/blob/main/scripts/figure4.py) The NIRISS-SOSS observations enable the detection of clouds in the atmosphere of WASP-39b.Clear atmosphere models cannot explain the amplitudes of all of the water vapour features simultaneously, which strongly indicates the presence of clouds (see Methods and Extended Data Fig. 6).The atmospheric models explored here indicate the presence of non-gray and non-homogeneous clouds, with model preferences of and greater for models with both 8σ non-gray and non-homogeneous clouds over models with gray homogeneous clouds only.This model preference is driven by the decrease in transit depth between 2-2.3 m (see Extended Data µ Fig. 7, panel a), which cannot be explained by gray clouds uniformly distributed along the terminator (see Extended Data Fig. 7, panel a).Moreover, within the various cloud treatments tested here (gray, gray + power-law, flux-balanced clouds, see Methods), both parametric and droplet sedimentation models indicate a preference for inhomogeneous cloud coverage of % around the planetary day-night terminator because it better explains the decrease ∼ 50 − 70 in transit depth between 2-2.3 m.These results confirm the necessity for broad wavelength µ coverage, and information regarding the cloud coverage of the planet, to constrain the atmospheric metallicity of a planet 2,4,11 Atmospheric circulation and cloud microphysical models have predicted that the cloud structure varies significantly along the terminators of hot Jupiters [26][27][28] .In particular, different compositions of clouds have different condensation temperatures and thus likely have different cloud coverage at the terminator 26 .Further studies combining temperature difference of east/west terminators to microphysical cloud models may be able to use the measured cloud coverage to determine the cloud composition of WASP-39b.Previous indications of non-gray or non-homogeneous clouds 5,[29][30][31][32][33] have relied on a single or small number of spectroscopic points, making our inference here for WASP-39b of non-gray cloud with inhomogeneous terminator coverage in the transmission spectrum of an exoplanetary atmosphere the most confident to date.These constraints on the physical properties of clouds, alongside the multiple spectral features across a broad wavelength coverage, are key to breaking well-known degeneracies between the metallicity and cloud-cover in atmospheric models 2,4,11 and deriving constraints on the bulk atmospheric properties.
The precision and wavelength coverage of NIRISS-SOSS should allow us to obtain more precise and robust constraints on atmospheric composition and tracers of planet formation than most previous transmission spectroscopy observations.The super-solar metallicity of WASP-39 b and the solar-to-super-solar K/O are in agreement with previous studies of mass-metallicity trends in transiting exoplanets 18,22,34,35 .If confirmed with further detailed modelling, a super-solar K/O ratio in WASP-39 b's atmosphere would likely indicate enrichment due to the accretion of planetesimals [20][21][22] , although the measurements of potassium and oxygen abundances for the host star are also needed in order to establish this result.Similarly, the suggestion of sub-solar C/O and super-solar metallicity may be compatible with a planetesimal accretion scenario e.g., 36,37,38 .The combination of a super-solar metallicity, super-solar K/O ratio, and sub-solar C/O ratio may suggest the planet formed beyond the H 2 O snowline followed by inward migration, for which theory predicts efficient accretion of planetesimals at -AU e.g., 20,39 .At those orbital ∼ 2 10 distances, the planetesimals likely contain K rock and H 2 O ice but almost no CO ice e.g., 40,41 , which explains the sub-solar C/O and super-solar K/O ratios along with a super-solar metallicity if a sufficient amount of planetesimals were accreted.However, fully understanding the possible formation pathways of this planet requires statistical constraints on the complete chemical inventory of the planet and the relative abundances of the carbon-oxygen-and alkali-bearing species.Such efforts will be possible when applying retrieval techniques to the complete transmission spectrum of WASP-39b from 0.5 to 5.5 m that is being produced by the Transiting µ Exoplanet Community ERS Program.Our results validate JWST's NIRISS-SOSS as an instrument mode capable of producing exquisite exoplanet atmosphere measurements and major advancements in our understanding of planet formation and exoplanet atmospheres.

Data Reduction
We applied six independent data reduction and light curve fitting routines to the data: nirHiss, supreme-SPOON, transitspectroscopy, NAMELESS, iraclis, and FIREFly.Below, we first describe the major reduction steps taken by each, followed by their light curve fitting methodologies.We note here that in each pipeline the position of the SOSS trace was found to match near-perfectly with that measured during commissioning.Each pipeline therefore used the default wavelength solution for SOSS measured during commissioning and provided in the jwst_niriss_spectrace_0023.fitsreference file.We present a summary of all pipelines in Extended Data Table 1.
The nirHiss Pipeline nirHiss is a Python open-source package that uses the Stage 2 outputs from the Eureka!pipeline and performs additional background and cosmic ray removal as well as extraction of the stellar spectra.For this analysis, we took the uncalibrated images and ran our own Stages 1 and 2 calibration using Eureka! 42, an open-source package which performs spectral extraction and light curve fitting for several JWST instruments.We use the default steps presented in Eureka!, which includes detector-level corrections, production of count rate images, application of physical corrections, and calibrations to individual exposures.
Next nirHiss removes background noise sources in a multi-step process.The zodiacal background is first removed by applying the background model provided on the STScI JDox User Documentation website. 1 The background is scaled to a small region of each science integration where there was no contamination from any of the orders; in this case, x [190,250], ∈ y [200,500].The average scaling-calculated here to be 0.881-is applied to all science ∈ integrations.Second, a model of 0 th order contaminants is built using the F277W integrations.The F277W integrations were taken after the transit of WASP-39b with the GR700XD/CLEAR pupil element and the F277W filter (throughput centred at m with a bandwidth of λ = 2. 776µ m).These observations consist of ten integrations with an exposure time of 49.4s.λ = 0. 715µ Observations with the F277W filter contain only the spectral trace of order 1 in the region where pixels, thus allowing for the detection and modelling of 0 th order contaminants across  ≤ 460 the majority of the detector.A median F277W frame is created to identify and mask any bad data quality pixels.
To ensure no additional noise is added from the F277W median frame, we create a 2D background model map using photutils.Background2D.To identify regions of the background, we masked the upper-left corner, where the trace is located, and any regions with , which includes the 0 th order sources.For photutils.Background2D, we used a filter > 1. 5σ size of (3,2) pixels and a box size of (2,2) pixels.Once the background is removed from the median F277W frame, we apply a Gaussian filter with a width of 2 to smooth out any additional small-scale background noise.To apply the median F277W frame to the Stage 2 science integrations, we scaled it to two isolated 0 th order sources in the science integrations at  1 ∈ [900,1100], [150, 250] and [1800,2000], [150, 250].We applied the average scaling to all integrations.We found the average F277W background scaling to be 2.81.We apply the scaled background frame to each time-series integration (TSO) integration.
Once the 0 th order contaminants are removed, we trace the location of Orders 1 and 2. The spatial profile for NIRISS-SOSS along the column is double-peaked, with a slight dip in the middle.We developed a routine to identify the trace locations using a three-step approach to identifying each order.For each column in the first order trace, we identify the locations of the two peaks, or "ears" and assume the middle of the trace is the median row-pixel between the two "ears".We repeat this process for the third and second orders in that sequence, masking orders once they have been traced.We chose to identify the third order before the second order since it is better spatially resolved and does not overlap with any other orders.The routine creates one main set of traces from a median frame of all observations which is used to extract the stellar spectra.As an additional output, we track the changes in (x,y) pixel positions of each order on the detector across all integrations.
After the traces are identified, we continue our reduction to remove any additional noise and cosmic rays/bad pixels.We perform additional 1/ noise correction following the routine  presented in transitspectroscopy (described below).Finally, nirHiss identifies and interpolates over cosmic rays.To identify cosmic rays we used the L.A. Cosmic technique wrapped into ccdproc 43,44 , which identifies pixels based on a variation of the Laplacian edge detection.We identify cosmic rays as pixels with via this method.We interpolate over any σ > 4 additional bad pixels by taking the median value of the two surrounding pixels along the column.We extract the spectra using a box extraction routine, and ignore any contaminants from overlapping orders or from any potential background orders.We use a box diameter of 24 pixels for both Orders 1 and 2. We do not extract the spectra from Order 3.

Extended Data Figure 1: Comparison of (x,y) position of NIRISS Orders 1 and 2 across the detector as modelled from different reduction pipelines. (a/c):
Each pipeline traces the curvature of Orders 1 and 2 using different methods.We show the best-fit trace for Order 1 in the top row, and Order 2 in the bottom row.We highlight zoomed-in regions to further examine differences.We note that iraclis uses the JWST-provided spectral trace.There is generally good agreement between the traces across the entire detector (<1 pixel deviations), with the strongest deviations towards the ends of each trace (e.g., x pixel position < 500 for Orders 1 and 2. (b/d): We provide an example spatial profile along the column for Order 1 (b) and Order 2 (d) at x = 1250.(https://github.com/afeinstein20/wasp39b_niriss_paper/blob/main/scripts/edfigure1.py)

The supreme-SPOON Pipeline
In parallel, we reduce the WASP-39 b TSO with the independent supreme-SPOON (supreme-Steps to Process sOss ObservatioNs) pipeline, which processes SOSS TSOs from the raw, uncalibrated detector images to extracted 1D light curves.An outline of the specific steps is presented below.
For detector-level processing supreme-SPOON closely follows Stage 1 of the jwst pipeline.All default steps, up to and including the reference pixel correction, are run using their default settings.The reference pixel step is known to provide an inadequate correction of 1/ noise for SOSS observations; however, we include it in order to remove group-to-group variations in the bias level, as well as even-odd row variations.At this stage, we remove the zodiacal background from each group.This is accomplished by first calculating a group-wise median frame, and scaling the model background provided in the STScI JDox to the flux level of each group in this median -yielding eight background models, one for each group.The region chosen to calculate the scaling was x [300,500], y [210,250], where there is minimal ∈ ∈ contamination from any of the SOSS orders.The n th background model is then subtracted from the corresponding group of each integration.
We then proceed to a more in-depth treatment of 1/ noise.Unlike the other pipelines  used in this work, supreme-SPOON treats 1/ noise at the group level instead of the integration  level.1/ noise is a time-varying noise source introduced by the voltage amplifiers during the  readout of the detector, and therefore the 1/ pattern will vary from group-to-group, even within  a given integration.To perform the 1/ correction, first a median out-of-transit frame is  calculated for each group.This group-wise median is then scaled to the flux level of each frame in a given group via the transit curve, and subtracted off -revealing the characteristic 1/ striping in the residuals.A column-wise median of this residual map is then subtracted from the original frame.The trace residuals as well as any bad pixels are masked in the median calculation.
From this point we once again proceed with the standard Stage 1 steps of the jwst pipeline, with the exception of the dark_current step, to obtain the supreme-SPOON Stage 1 outputs.The dark current subtraction step is skipped as it was found to re-introduce 1/ noise  into the data.The dark current level is additionally extremely small (several 10s of electrons/s compared to many thousands for the target signal) and can thus be safely ignored.supreme-SPOON only applies the assign_wcs, srctype, and flat_field steps of the Stage 2 jwst pipeline to the Stage 1 products.The background subtraction was already performed as part of Stage 1 calibrations.Furthermore, the flux calibration steps (pathloss which accounts for light incident on the telescope primary mirror which falls outside of the SUBSTRIP256 subarray, and photom which performs the actual photometric flux calibration) are skipped, both because an absolute flux calibration is unnecessary for relative spectrophotometric measurements, and a wavelength-dependent flux calibration is nonsensical for SOSS where contributions from multiple wavelengths from all orders impact a single pixel.At this point, supreme-SPOON identifies any remaining hot pixels via median filtering of a median stack of all frames and interpolates them via the median of a surrounding box.These products are the supreme-SPOON Stage 2 results.
Stage 3 of the supreme-SPOON pipeline is the 1D extraction.This can be performed via two different methods: the first is a simple box aperture extraction on each order, ignoring the order contamination.The second uses ATOCA (Algorithm to Treat Order ContAmination) 45 to explicitly model the order contamination.Briefly, ATOCA constructs a linear model for each pixel on the detector, including contributions from the first and second diffraction orders, allowing for the decontamination of the SOSS detector -that is, ATOCA constructs models of both the first and second orders individually, thereby allowing a box extraction to be performed on each free from the effects of order contamination.Although the effects of this order contamination for differential measurements (e.g., exoplanet atmosphere observations) are predicted to be small (∼ 1% of the amplitude of the expected spectral features) 45,46 , in the quest to obtain the most accurate possible transmission spectra, this contamination effect is important to take into account.ATOCA is currently built into the extract1dstep of the official jwst pipeline, although it not currently the default option and must be toggled to on via the "soss_atoca" parameter.To improve the performance of ATOCA, we do not use the default specprofile reference file included in the jwst pipeline, but instead construct estimates of the underlying spatial profiles of the first and second order, upon which ATOCA relies, using the APPLESOSS (A Producer of ProfiLEs for SOSS) algorithm 46 .We determine the centroid positions for each order on a median stack using the "edgetrigger" algorithm 46 , and these positions are found to match to within a pixel with the default centroids contained in the jwst_niriss_spectrace_0023.fitsreference file; the spectrace file is available on the JWST Calibration Reference Data System (CRDS). 2 The SOSS trace position is furthermore highly stable over the course of this TSO, with RMS variations in x and y positions of 5 mpix, and RMS rotation of 0. The transispectroscopy Pipeline This third pipeline analysis combines the jwst pipeline Stage 1 "rateints.fits"files with transitspectroscopy 47 .transitspectroscopy completes stellar spectral extraction as well as transit fitting.
The trace positions for NIRISS orders 1 and 2 were determined using transitspectroscopy.trace_spectrum.This routine cross-correlates an input function with each column in the detector to find the centre of the different traces via the maximum of the resulting cross-correlation function.In order to follow the shape of the NIRISS order profiles, an input function consisting of a double Gaussian was used with parameters that were trained on the NIRISS/SOSS observations of HAT-P-14 b (JWST Program ID 1541; PI Espinoza): . The trace for order 2 was not fit for pixels , µ 2 =− 7. 5; σ 1 = 3. 0; µ 2 = 7. 5; σ 2 = 3. 0 ≤ 1040 as the throughput is not high enough for the method to robustly fit the trace without incorporating nearby contaminants.After identifying the trace positions with this method for both Order 1 and Order 2, both traces are smoothed using a series of spline functions.We find the best-fit parameters for Order 1, which were trained on the HAT-P-14 b observations are: The zodiacal background was removed by scaling the model background provided on the STScI JDox User Documentation.This model was compared to a small region of the median science integrations where there was little to no contamination from the orders (x [500,800], y ∈ ∈ [210,250]).The ratio of all the pixels in this region versus the pixels in the background model was computed, ordered, and the median ratio of all the second quartile pixels was used as the scaling factor between the background model and the data, which was found to be .All the 0. 909 integrations had this scaled background subtracted.
Each integration is corrected for 1/ noise with the following procedure.First, all the  out-of-transit, background-corrected integrations are median combined and scaled by the relative flux decrease produced by the transit event at each integration (i.e., 1.0 for out-of-transit integrations, or about 0.976 for mid-transit).These scaled median frames are then subtracted from each individual integration, which then leaves in the frame only detector-level effects such as 1/ noise.We then go column by column and take the median of all pixels in this residual  frames within a distance of 20 to 35 pixels from the centre of the trace, and use this as an estimate of the contribution from 1/ noise to that given column.This value is then removed  from each pixel within 20 pixels from the trace on that column.No correction for Order 1 contamination on Order 2 was made as the contribution is negligibly small in this case 45 similarly for Order 1 contamination in Order 2 in our extraction.
To extract the resulting background and 1/ -corrected spectrum, we used the  transitspectroscopy.spectroscopy.getSimpleSpectrumroutine with a 30-pixel total aperture for both orders.In order to handle obvious outliers in the resulting spectrum due to, e.g., uncorrected cosmic rays and/or deviating pixels, we median-normalised the spectra for each integration and combined them all to form a "master" 1D spectrum for both orders 1 and 2. The median was taken at each wavelength, as well as the error on that median, and this was then used to search for outliers on each individual integration at each wavelength.If outliers were 5σ found, they were replaced by the re-scaled version of this median "master" spectrum.
Extended Data Figure 3: Example extracted stellar spectra from different reduction pipelines.The inset axes highlight the peak of the spectra.The supreme-SPOON spectra are scaled by a factor of 72 to compare the overall shape of the spectra, rather than the extracted flux counts.Top: The extracted stellar spectra from Order 1.All pipelines are in relatively good agreement, while the shape of the iraclis data changes slightly at λ < 1.1μm.This is likely due to different traces which were used in the spectral extraction routine.Bottom: The extracted stellar spectra from Order 2. Differences at λ = 0.725 and 0.86μm are due to differences in removing zeroth order contaminants in the background.The \texttt{iraclis} pipeline does not extract data out past λ > 0.85μm, which is where the order overlap region begins.Across all pipelines, the shape of the spectra, as well as overall absorption features, cosmic ray removal techniques, and noise levels are consistent.(https://github.com/afeinstein20/wasp39b_niriss_paper/blob/main/scripts/edfigure3.py) The NAMELESS Pipeline Starting from the jwst pipeline Stage 1 products, we use the NAMELESS (Niriss dAta reduction MEthod for exopLanEt SpectroScopy) pipeline to go through the jwst pipeline Stage 2 with the addition of custom correction routines.
First, we go through the assign_wcs, srctype, and flat_field steps of the jwst pipeline Stage 2, opting for a custom background subtraction routine and skipping the pathloss and photom steps as absolute flux calibration is not needed.After flat-field correction, we scale the model background provided on the STScI JDox User Documentation to a region of the median frame where the contribution from the tail of the three orders is lowest (x [200,250], y [400,600]).From the distribution of the scaling values of all pixels within the ∈ ∈ defined region, we take the 16 th percentile as our scaling value and subtract the scaled background frame from all integrations.
We subsequently correct for 1/ noise by performing a column-by-column subtraction for  each median frame subtracted integrations.The median frame is computed from the out-of-eclipse integrations (integration # [200,400]) and scaled to each individual integration ∈ by dividing the sum of the pixels within the first order by that of the median frame.We then subtract the scaled median frame from all integrations, perform the column-by-column subtraction on the residual frames, and add back the scaled median frame to the corrected residual frames to obtain the 1/ corrected integrations.

𝑓
We detect outliers frame-by-frame using the product of the second derivatives in the column and row directions.This method works particularly well for isolated outliers, as this leads to a strong inflection that corresponds to a large second derivative.Because the spectral orders also lead to larger second derivative values, we divide the frames into windows of 4 4 × pixels, compute the local second derivative median and standard deviation, and flag any pixel that is more than four standard deviations away from the median.Furthermore, we also flag pixels with null or negative flux.All identified outliers are set equal to the median value of the window it was identified in.
Finally, we proceed with spectral extraction of the corrected frames by first tracing the sections of the spectral orders that we wish to extract.We trace orders 1 and 2 from x 1 [4,2043]  ∈ and x 2 [4,1830] respectively.The centre of the traces is found for each individual column by ∈ performing a convolution of the profile with a gaussian filter, where we use the maximum of the convolved profile as the centre of the trace.For the tracing of the second order, we keep the centre of the trace fixed below x as the flux from the first order can bias the tracing = 900 method.Furthermore, we smooth the positions of the trace centroids using a spline function with 11 and 7 knots for the first and second order respectively.We perform spectral extraction of the first and second orders at all integrations using the transitspectroscopy.spectroscopy.getSimpleSpectrumroutine with an aperture width of 30 pixels.

The iraclis Pipeline
We used the jwst pipeline Stage 1 "rateints.fits"files with modified routines from iraclis 48,49 , which was initially designed for HST.The modified routines will be part of the version 2 of iraclis, which will become publicly available in the near future.The routines applied to the "rateints.fits"files were flat fielding, bad pixels and cosmic rays correction, sky background subtraction, 1/ noise correction, X-and Y-drifts detection, light curve extraction,  light curve modelling and planetary spectrum decontamination.
We started our analysis by dividing the images by the appropriate flat field frame (jwst_niriss_flat_0275.fits), as provided by the JWST CRDS.The next step was the bad pixels and cosmic rays correction.As bad pixels we used those with a positive DQ flag in the "rateints.fits"files, excluding the warm pixels, as their large number did not allow for a reliable correction.We also identified extra outliers (cosmic rays or other artefacts) by calculating two flags for each pixel: the difference from the average of the ten horizontally neighbouring pixels (x-flag) and the difference from the average of the ten vertically neighbouring pixels (y-flag).If a pixel's x-flag was 5 larger than the other pixels in the column and its y-flag 5 larger than the σ σ other pixels in the row, it was identified as a cosmic ray (see also 49 ).Both bad pixels and outliers were replaced with the value of a 2D interpolation function, created from the rest of the pixels, similarly to analyses with HST 49 .
We then subtracted a column-based sky background frame and a column-based 1/ noise  frame from each integration.For each image, we first used a trace filter (value 0.001 in the > jwst_niriss_spectrace_0023.fits,provided by the JWST CRDS), and a column-based 1 median × absolute deviation filter to find the illuminated pixels.Then, we calculated the column-based median of the image-using only the unilluminated pixels-and subtracted it from the image.Finally, we calculated the column-based median of the IMFD (Image-MedianFrame Difference)-using only the unilluminated pixels-and subtracted it from the image.This process is not efficient in subtracting 100% of the background contamination, which was removed during the last analysis step (spectrum decontamination).
After reduction, X-and Y-drifts were detected relatively to the first image by comparing the sums along the columns and the rows, respectively, similarly to HST 49 .We then extracted the light curves at high resolution, to avoid the 1/ noise.For each spectroscopic image, we initially  divided each pixel into a 100 100 grid of subpixels, and, for each subpixel, we calculated the × distance from the trace (CD) and the wavelength ( ), creating the and the , λ   λ  respectively.was assigned to each subpixel directly from the wavelength solution (interpolated λ wavelength solution from the jwst_niriss_wavemap_0013.fits file, provided by the JWST CRDS, shifted by the detected X-and Y-drifts).CD was calculated as the distance between the centre of the subpixel and the point of the trace with the same (interpolated trace from the λ jwst_niriss_spectrace_0023.fitsfile, provided by the JWST CRDS, shifted by the detected Xand Y-drifts).Our high-resolution bins had a -width of 10 Å, ranging between 0.62 and 0.85 m λ µ for order 2 and between 0.85 and 2.8 m for order 1, and a -width of 1.5 pixels, ranging from µ  -25 to 25 pixels.
Finally, to construct the light curve of each bin, we applied the following smoothed aperture mask on each spectroscopic image and summed the values of all the subpixels.We chose a smoothed aperture, similarly to HST to reduce the effects of jitter noise: where , , and are the bin boundaries and the smoothing factor along the  1  2 σ  cross-dispersion axis, , , and are the bin boundaries and the smoothing factor along the λ 1 λ 2 σ λ dispersion axis.For the smoothing factors we used the values of pixels and σ  = 0. 015 σ λ = 1 Å -i.e., 10% of the bin size.We chose these values for the smoothing factors because lower ∼ values would effectively create a sharp-edge aperture, while larger values would force the bins to overlap substantially.

FIREFly
While FIREFLy (Fast InfraRed Exoplanet Fitting for Light curves) 50 was written and optimised for reducing NIRSpec-PRISM and G395H time series observations, it worked well on the NIRSISS-SOSS dataset, where it selected and processed the spectrophotometry from the first order only with minimal tuning or intervention.FIREFLy is not written in such a way to extract data from Order 2 ( m).In our reduction, we perform standard calibrations on the raw λ < 0. 9µ data using the jwst Python pipeline, and we perform bad pixel and cosmic ray cleaning on each integration.We perform 1/ destriping and background subtraction using a pixel mask generated  from the temporally-medianed image that selects regions of the data image below a specified count threshold.We extract the spectrophotometry using an optimised aperture extraction width that is constant in wavelength.The aperture width is selected such that the scatter of the resulting out-of-transit white-light photometry is minimised.We bin the cleaned spectrophotometry in wavelength to create 120 variable-width spectral channels with roughly equal counts in each.

Light Curve Fitting and Transmission Spectra
We used a suite of light curve fitting routines to fit the extracted light curves.Each routine fit for orbital parameters from the broadband white-light curves per each order (see Extended Data Figure 4).For the spectroscopic light curves, most routines (nirHiss/chromatic_fitting, supreme-SPOON/juliet, transispectroscopy/juliet, and NAMELESS/ExoTEP) fixed the orbital parameters (that is the mid-transit time, , semi-major axis to stellar radius ratio  0 , impact parameter , eccentricity ) to the same values to ensure consistency.These / *   parameters were fixed to their best-fitting values from the transitspectroscopy/juliet white light curve fit, except for which was fixed to the value obtained from the white light  0 curve in each case.The orbital period of WASP-39 b was fixed to 4.05528 d 9 for all fits.This left the planet-to-star radius ratio , the limb-darkening coefficients, and parameters for any   / * additional systematics models to vary.These four routines also fit spectroscopic light curves at the native instrument resolution.However, two routines, iraclis and FIREFly instead fixed the orbital parameters in their spectroscopic fits to values obtained via their white light curve fits.iraclis also fit directly for the orbital inclination, , as opposed to and like the other   / * routines.Additionally, iraclis and FIREFly fit for their spectrophotometric light curves at the pixel resolution.We present all of the best-fit white-light curve parameters for Order 1 in Extended Data Table 2. Additionally, for the spectroscopic light curve fits, we only considered the region of order 2 with wavelength 0.85 µm, as the 0.85 -1.0 µm range is covered at higher < signal-to-noise by order 1.

chromatic_fitting
chromatic_fitting is an open-source Python tool for modelling multi-wavelength photometric light curves.This tool is built on the framework of chromatic, a package for importing, visualising and comparing spectroscopic datasets from a variety of sources, including Eureka! and the JWST pipeline.In this paper we applied chromatic_fitting to the nirHiss reduction.
chromatic_fitting utilises the PyMC3 (NUTS) sampler 51 to fit the exoplanet transit model to the light curves.First, we fit the white light curves for order 1.The white light curve was generated using an inverse variance weighted average of the unbinned data.We fixed the orbital period to 4.05528 d 9 and assumed a circular orbit.We fit for the mid-transit epoch , the  0 stellar mass and radius , the impact parameter , the planet-star radius ratio implementation we ran 4000 tuning steps and 4000 draws for the white light curve and the mean parameter values are shown in Extended Data Table 2.We checked for convergence using the rank normalised R-hat diagnostic 55,56 .
For each spectroscopic light curve we fixed the period , transit epoch , eccentricity ,   0  semi-major axis in stellar radii and impact parameter to the white-light parameter values / *  from the transitspectroscopy/juliet routine (Extended Data Table 2).We then fit for the planet-star radius ratio , quadratic limb-darkening coefficients ( , ) and out-of-transit -for all of these parameters we assumed the same normal priors as for the white light  0 curve.We also included a second-order polynomial in time with the same priors as the white light curve fit.For each wavelength we ran 2000 tuning steps and 2000 draws.The final transmission spectrum was taken as the mean value drawn from the posterior distribution for the planet-star radius ratio with 1 uncertainties extracted from the 16-84 th HDI (highest density σ interval) region.

juliet
We applied the juliet package 57 for light curve fitting to the products of multiple reduction pipelines described above.Here, we give a general overview of the methods and include exact details per fit when appropriate.
For the supreme-SPOON reduced stellar spectra, we fit the white light curve for the mid-transit time, , the impact parameter, , the scaled orbital semi-major axis, , and the  0  / * scaled planetary radius, ; assuming a circular orbit.We also fit two parameters of a   / * quadratic limb-darkening model following the parameterization of 58 , as well as an additive scalar jitter and the two parameters of a linear trend with time.We therefore fit seven parameters to the white light curve for each order, using wide, flat priors for each case.We then proceeded to fit the light curves from each individual wavelength bin at the native detector resolution -that is two pixels per bin to roughly account for the extent of the PSF in the spectral direction.This results in 1020 bins for Order 1, and 520 bins for Order 2 as we only consider wavelengths < 0.85 µm.For the spectroscopic fits, we fixed the central transit time to the white light curve value, and the other orbital parameters as described for chromatic_fitting.For the linear trend with time, we used the white light posterior for each of the two parameters as prior distributions for all wavelength bins, whereas for the limb-darkening parameters, we adopted a Gaussian prior centred around the predictions of the ExoTiC-LD package 54 with a width of 0.1.As the SOSS throughput files included with ExoTiC-LD did not cover the full wavelength range of both orders, we instead used the throughputs determined during commissioning and included in the spectrace reference file of the jwst pipeline.We truncated the Gaussian prior at 0 and 1, to prevent the limb-darkening parameters from straying into unphysical regions of the parameter space.We then used flat, uninformative priors for the remaining two parameters, the scaled planetary radius and the scalar jitter.The supreme-SPOON white light curve fits have for Order 1, and for Order 2.
For the transitspectroscopy reduced stellar spectra, we first fit the white-light light curves of Order 1 and 2 separately.For these, as suggested above, the period was fixed, but all the rest of the parameters were allowed to vary.In particular, we set a normal prior on the time-of-transit centre of days, where the first value denotes the mean and  2459787.5, 0. 2

( )
the second the variance of the prior.A normal prior was also set on and / * ∼  11. 37, 0. 5

( )
where the means were set following the work of 59 , but the variances are large to account for the variation of these parameters in the literature between different authors.We set a uniform prior for the planet-to-star radius ratio between and , and fixed eccentricity to 0. In addition to 0 0. 2 those, we fit for a mean out-of-transit offset with a normal prior of , and a jitter term  0, 0. 1

( )
added in quadrature to the error bars with a log-uniform prior between 10 and 1000 ppm.To account for systematic trends in the data, we use a Gaussian Process (GP) via celerite 60 with a simple Matèrn 3/2 kernel to parametrize those trends.We set log-uniform priors for both the amplitude of this GP from to 1000 ppm, and for the time-scale from days to days.10 We use the framework of 58 to parameterize limb-darkening via a square-root law, which, following 61 , is one of the laws that should give the best results at this level of precision.
For the wavelength-dependent light curves we used a similar setup with two main differences.The first is that we fix the time-of-transit centre, , and to their white light / *  values.The second is that we set truncated normal priors on the transformed limb-darkening coefficients between 0 and 1, with standard deviations of and means obtained via  1 ,  2 ( ) 0. 1 the following method.First, we obtain the non-linear limb-darkening coefficients using an ATLAS stellar model with the closest properties to WASP-39's via the limb-darkening software 62 .Then, the square-root law limb-darkening coefficients are obtained following the algorithm of 63 , and those are transformed to the parameterization using the equations  1 ,  2 ( ) in 58 ; these are then set as the mean for each wavelength-dependent light curve.We note that we fit the light curves at the pixel level, which means fitting one light curve per detector column.We fit them in parallel using the transitspectroscopy.transitfitting.fit_lightcurvesroutine.

ExoTEP
For the NAMELESS reduction, we perform light curve fitting on the extracted spectrophotometric observations using the ExoTEP framework 64 .We first fit the white light curves of both orders 1 and 2 separately.We fit for the mid-transit time , the planet-to-star  0 radius ratio , and quadratic limb-darkening coefficients (u , u ) 65,66 , while fixing the impact parameter and semi-major axis to the values of the best order 1 white light curve fit from  / * the transitspectroscopy/juliet analysis.We also fit for the scatter , as well as a linear σ systematics model with an offset and slope .Uniform priors are considered for all parameters.
We additionally only discard the first 10 minutes of observations (10 integrations) to remove the exponential ramp.For all light curves, we compute the rolling median for a window size of 11 integrations and bring any data point that is more than four standard deviations away from it to the median value.We fit the light curves using the Markov chain Monte Carlo Ensemble sampler emcee 67 for 1,000 steps using four walkers per free parameter.The first 600 steps, 60% of the total amount, are discarded as burn-in.We then fit the spectroscopic light curves, keeping  0 fixed to its white light value, at a resolution of three pixels per bin for order 1 (680 bins) and one pixel per bin for order 2 from 0.6-0.9m (675 bins).We used 1,000 steps for the spectroscopic ∼ µ fits, once again discarding the first 600 as burn-in.

iraclis
We analysed all the light curves using the open source Python package PyLightcurve 68 .For every light curve, PyLightcurve: (I) calculates the limb-darkening coefficients using the ExoTETHyS package 69 , the wavelength range of the bin, the response curves for each of the NIRISS orders (jwst_niriss_spectrace_0023.fitsfile, provided by the JWST CRDS), and the stellar parameters ( =5540 K, =4.42 cm/s , Fe/H=0.14 dex 70 ); (II) finds the    2 maximum-likelihood model for the data (an exposure-integrated transit model together with a quadratic trend model using the Nelder-Mead minimisation algorithm included in the SciPy package 71 ; (III) removes outliers that deviate from the maximum-likelihood model by more than three times the standard deviation of the normalised residuals; (IV) scales the uncertainties by the RMS of the normalised residuals, to take into account any extra scatter; (V) and, finally, performs an MCMC optimisation process using the emcee package 67 .We initially modelled the first order white light curve (sum of all bins above 0.85 m with out-of-transit fluxes above 20 µ DN s ) and fit for the white , the orbital parameters, and , and the transit mid-time.uniform contamination.The uniform contamination fixes issues of sky background over-or under-correction.It also corrects for order overlap.After the decontamination described above, there was still a contaminating source affecting the spectrum around 0.72 m, which was not µ uniform due to the PSF.To separate this source we applied Independent Component Analysis (ICA) on the stellar spectra extracted from various distances from the trace.We used two components to describe the contaminating source and one to describe the stellar spectrum.
Finally, we estimated the per wavelength bin using the weighted average of all the bins   / * ( ) 2 that had the same wavelength range.We only took into account the bins that had out-of-transit fluxes above 20 DN s .This choice effectively applied an optimal aperture size for each −1 wavelength bin.

FIREFly
To extract the transmission spectrum, we fit for the planet's transit depth at each wavelength using a joint light curve and systematics model.We use the orbital parameters recovered from an MCMC fit to the white light curve, and fix them at each wavelength channel for our fit.We fit for the two quadratic limb-darkening terms and at each wavelength   channel.We find that the best-fit limb-darkening coefficients are uniquely determined, and deviate by a constant offset relative to model coefficients.Our fits are performed iteratively using the Python package lmfit.The light curves show a typical photometric scatter of 0.3% per integration, and the typical transit depth uncertainties vary between 150-300 ppm, which is in line with near photon-limited precision.More details of the FIREFly fitting routine can be found in 50 .
Extended Data Figure 5: Transmission Spectra for WASP-39b for all reduction techniques.Our best-fit reference model to the nirHiss spectrum (red) is plotted as a black solid line on both panels and the spectra are separated into three panels for ease of reading.Wavelengths which overlap with 0 th order contaminants are masked.(https://github.com/afeinstein20/wasp39b_niriss_paper/blob/main/scripts/edfigure5.py)

Atmospheric Models
To interpret the measured transmission spectrum, we performed an extensive comparison with grids of synthetic transmission spectra.We tested several independent atmospheric models to avoid any model-dependent interpretation of the data.Unless otherwise noted, all of our grids have assumed radiative-convective-thermochemical equilibrium to estimate atmospheric compositions.The exploration of atmospheric models with fewer assumptions (e.g., without the assumption of chemical equilibrium with metallicity and C/O as free parameters) and those considering other effects of disequilibrium chemistry is left for future work.
We derive basic interpretations for the observed spectrum based on four independent model grids, ATMO, PHOENIX, PICASO, ScCHIMERA .Each grid contains pre-computed transmission spectra at various atmospheric properties, such as metallicity [M/H], carbon-to-oxygen C/O ratio, and cloud properties, using from gray to Mie-scattering cloud opacity (see next subsection for details).The ScCHIMERA grid considers additional model advancements 1) various cloud treatments, including a gray cloud, gray + power-law cloud opacity, and physically motivated (i.e., droplet sedimentation) cloud model, 2) the impact of inhomogeneous cloud coverage along the planetary terminator, and 3) the K/O ratio as a grid dimension.ScCHIMERA provides the best fit to the observations compared to the other three grids and informs the results presented in the main text.

Grid Search with Pre-Computed Forward Models
Here, we introduce the independent grids of pre-computed transmission spectra, their model description, and the main results from these grid fits.We first present the three grids that do assume horizontally homogeneous clouds.

PHOENIX
The atmospheric PT profile is computed using the 1D radiative-convective equilibrium model PHOENIX [81][82][83] 84 , and atomic lines from the database of Kurucz 85 .For cloudy models, the small non-gray cloud particle opacity is treated as a sum of Rayleigh scattering opacity of all gas species enhanced by a factor of either (clear atmosphere) or ; 0 10 large gray particle opacity is treated as gray cloud deck pressure levels of , , and 10 mbar.In 0. 3 3 total, the PHOENIX grid consists of 95 cloud-free and 380 cloudy atmosphere models.We only consider horizontally homogeneous clouds in the PHOENIX grid fits.

PICASO 3.0
Similarly to the grids of models presented above, we precomputed atmospheric pressure-temperature (PT) profiles using the 1D radiative-convective equilibrium model PICASO 3.0 [86][87][88][89]  , CO 93 , K from VALD3 80 .For cloudy models, we post-processed the computed PT profiles using the droplet sedimentation model Virga 94,95 that determines the vertical distributions of cloud mass mixing ratio and mean particle size from the balance between downward mass flux of gravitational sedimentation and upward mass flux of eddy diffusion.We vary vertically-constant eddy diffusion coefficients of K zz = 10 5 , 10 7 , 10 9 , 10 11 and vertically-constant sedimentation parameters of f sed =0.6, , , , .The f sed value is 1. 0 3. 0 6. 0 10. 0 defined as the ratio of the mass-averaged sedimentation velocity of cloud particles to the mean upward velocity of the atmosphere, and a smaller f sed yielding more vertically extended clouds 95 see e.g., 96 .We have assumed horizontally homogeneous clouds and accounted for the formation of MgSiO , MnS, and Na S clouds.Then, we postprocessed the atmospheric properties to compute synthetic transmission spectra.We note that the optical properties of the flux-balanced cloud is computed by the Mie theory 97 under the assumption of a log-normal particle size distribution with a mean particle size translated from f sed

94
. In total, the PICASO grid consists of 192 cloud-free and 3840 cloudy atmosphere models.
We compare the NIRISS SOSS spectrum (binned to R=300) to each of these model grids and summarise the best fits in the top panel of Extended Data Fig. 6.For each cloudy and clear model we tested, we compute between the data and the models, with χ (https://github.com/afeinstein20/wasp39b_niriss_paper/blob/main/scripts/edfigure6.py)

Grid Search with ScCHIMERA
The NIRISS transmission spectrum offers key insights into the atmospheric properties of WASP-39b over a broad wavelength range.The simultaneous detection of H 2 O and K, alongside possible indications of carbon-bearing species, allows us to explore equilibrium models for which the potassium-to-oxygen (K/O) ratio is an additional dimension besides the commonly employed C/O and metallicity parameters.Furthermore, as explained in the previous subsection (see also Fig. 3 demonstrating how clouds contribute to the NIRISS spectrum), the broad wavelength coverage of these NIRISS observations makes it possible to explore more complex cloud models beyond traditional gray and homogeneous cloud models.To explore these considerations, we implement the ScCHIMERA grid as explained below.

ScCIMERA
Previous implementations of this framework include [98][99][100] , where the methods are described in detail.Implementations of this procedure to JWST data include ref. 101 For a given set of planetary parameters, our methods pre-compute the temperature-pressure structure of the planetary atmosphere and the thermochemical equilibrium gas mixing ratio profiles.The computations are performed on a grid of atmospheric metallicity ([M/H], e.g., log 10 enrichment relative to solar 17 ) spaced at 0.125 dex values between 0 and 2.25 (e.g., 1 to 177 solar) and × carbon-to-oxygen (C/O) ratios at values of 0.2, 0.35, 0.45, 0.55, 0.65, 0.7, and 0.8.Unlike previous implementations of this framework, and to better understand the NIRISS-SOSS observations presented, we include a dimension to our grid exploring the potassium-to-oxygen ratio ([K/O], i.e., enrichment relative to solar 17 ) with spacing of 0.5 dex between -1 and 0,  10 and 0.1 dex between 0 and 1, overall spanning a range from -1 to 1, or 0.1 to 10 solar.In these × calculations, the atmospheric metallicity (M/H) scales the sum of K, C and O.This sum determines the final elemental abundances after scaling M/H, C/O and K/O.That is, the total oxygen elemental abundance is , the total carbon elemental abundance is ' = / /+/+1 , and the total potassium elemental abundance is .Additionally, ' = ' × / ' = ' × / we explore the energy redistribution (f) between the day and night sides of the planet 102 , with values of 0.657, 0.721, 0.791, 0.865, 1.0 , 1.03 , 1.12 , 1.217, 1.319 in our grid, where  = 1.0 and correspond to full day-to-night heat redistribution and dayside only redistribution, 2. 0 respectively.
The transmission spectrum of the planet is computed with CHIMERA 96,103,104 using the converged atmospheric structures.We compare the observations to these models in a Bayesian inference framework using the nested sampling algorithm MultiNest 105 through its python implementation PyMultiNest 106 , and obtain an optimal [M/H], C/O, [K/O], and f via nearest neighbour search in the grid., Na 111,112 , and K 111,113 , and were computed following the methods described in 114,115 .The cloud models considered are 1) a basic cloud model with a gray, uniformly vertically distributed cloud opacity ( ); 2) a gray+power-law κ  cloud model that accounts for non-gray opacity of small-size particles as a vertically uniform power-law opacity (i.e., a parameter for the scattering slope and and a Rayleigh-enhancement factor, e.g, 30,33,116,117 ) in addition to gray cloud component, which is expressed by gray cloud-deck of infinite opacity at a given atmospheric pressure; and 3) a droplet sedimentation model 94 (assuming enstatite grains) where parameters capture the eddy diffusion coefficient and the ratio of sedimentation velocity to characteristic vertical mixing velocity (see also the description of PICASO above).For cloud treatments 2 and 3, we also consider the possibility of inhomogeneities around the planetary limb by considering a linear combination of clear and cloudy models e.g. 11, key for breaking degeneracies between metallicity and cloud properties (e.g., 2,4 ).We assume the same PT profile for both cloudy and clear limbs in the inhomogeneous cloud models and leave an investigation on the possibility of different PT profiles on those regions to future studies.

Identification of Absorbers and Model Selection
We perform our Bayesian inference using all model combinations with the ScCHIMERA grid on four different data resolutions for the nirHiss transmission spectrum: R=100, R=300, native instrument resolution (R order 1 =910; R order 2 =830), and pixel-level resolution (R order 1 =1820; R order 2 =1660).Resolutions are given at the reference wavelengths of m for Order 1 and λ = 1.791µ m for Order 2. We test the robustness of our inferences against different binning and 0. 744µ convolution strategies and find the results, i.e., the bulk atmospheric properties [M/H], C/O, and K/O, consistent regardless of the resolution of the data.We find a fiducial combination of parameters that can best explain the spectrum (that we call the reference model) with full redistribution (f=1, matching predictions that planets in this temperature regime are unlikely to possess strong day-to-night temperature contrast [118][119][120] ), [M/H]=1.375(i.e., solar), ∼ 20 × C/O=0.2, and [K/O]=0.1.With these atmospheric properties, the data are best explained by the droplet sedimentation model (ScCHIMERA cloud model 3) and inhomogeneous cover.However, the gray+power-law model (ScCHIMERA cloud model 2) with inhomogeneous cover provides a comparable fit to the data.Using R=300 data, the homogeneous droplet sedimentation model (model 3) is preferred over homogeneous gray cloud (model 1) at , which strongly indicates ≳8σ the non-gray nature of cloud opacity.Meanwhile, the inhomogeneous droplet sedimentation model is preferred over the homogeneous droplet sedimentation model cloud at .This is 5σ evidence that for the same model 3 inhomogeneous cloud coverage is preferred.
We explore the contribution of different chemical species to our reference model by performing the Bayesian inference using the inhomogeneous cloud model 3 and artificially disabling the contribution of a selected chemical species, one at a time.By redoing the Bayesian inference, we are able to compare the Bayesian evidence by computing the Bayes factor and converting to a 'sigma' detection significance using the prescription in 2 (see also 22 ).We detected H 2 O at , > 30σ K at , CO at , and no significant detections of Na, CH 4 , CO 2 , HCN, and H 2 S. The 6. 8σ 3. 6σ best-fit metallicity across all models is -solar, the best fit K/O ratio -solar, and ∼ 10 30 × 1 2 × C/O ratio 0.2.Taking the average and standard deviation of the best-fit results for all 20 runs (i.e., 5 models on 4 data resolutions) we find an average M/H=19 solar with a standard × deviation of 5 solar and an average K/O=1.5 solar with a standard deviation of solar.× × 0. 26 ×

Wavelength Sensitivity to Inferences
We investigate the dependence of the inferred atmospheric properties on the spectral range of the observations by performing the same Bayesian inferences described above on the spectrum cloud models when using data blue-wards of m only.Without the information contained in 2 µ the dip in transit depth at 2.1 m, all cloud treatments provide an equally good fit and µ overestimate the transit depth between 2.0 m and 2.3 m. µ µ

K/O Inferences
We explore the possibility of constraining the potassium-to-oxygen ratio using NIRISS-SOSS.
As explained above, across different models and data resolutions, our results suggest that the observations of WASP-39b are best explained by a solar-to-super-solar K/O ratio.To further explore this, we repeat our Bayesian inference for all 20 model-data configurations (5 models × 4 resolutions) using the observations blue-wards of 0.8 m.From high-resolution to µ low-resolution observations and for all cloud model configurations, we find that all 20 runs prefer models with solar or super-solar K/O ratios for WASP-39b ranging from 1 to 10 solar.× The average across the 20 runs is 2.12 solar and a standard deviation of solar, with the × 2. 33 × relatively larger standard deviation resulting from two inferences of highly super solar K/O ratios of 7 solar or greater for observations at the pixel-level resolution.

Figure 1 :
Figure 1: Selection of spectrophotometric light curves and residuals for WASP-39b's transit observed with NIRISS-SOSS for Orders 1 and 2.An exoplanet transit model (solid line) was fitted to each light curve with chromatic_fitting using a quadratic limb-darkening law.The limb-darkening coefficients, planet-to-star radius ratio (R p /R * ), and out-of-transit flux were varied in each wavelength channel, while all other parameters were fixed.The residuals to the best-fit models are shown for each light curve.The wavelength range per each channel is denoted on the left plot, while parts-per-million (ppm) scatter in the residuals is denoted on the right plot.We calculate the ppm as the standard deviation of the out-of-transit residuals.The reductions are from the nirHiss and chromatic_fitting routines described in the Methods.(https://github.com/afeinstein20/wasp39b_niriss_paper/blob/main/scripts/figure1.py)

Figure 3 :
Figure 3: Interpretation of the constituents of the NIRISS WASP-39b transmission spectrum.(a/b): The top panel shows the comparison of WASP-39b's transmission spectrum from the nirHiss reduction (grey points) with respect to the best-fit reference model (black line).Each colored line removes a key constituent found in our best-fit reference model to demonstrate how the spectrum would change were these features not included.The removal of clouds, and H 2 O absorption from the reference model result in large-scale changes to the shape and depth of the transmission spectrum.Other sources of opacity with an impact on the spectrum are K, CO, and CO 2 .Residuals between the data and the reference model are plotted below.(c/d):The two bottom panels show the molecular absorption cross-sections for a selection of gases observable within the NIRISS bandpass.The middle panel highlights gases inferred by our analysis of WASP-39b's spectrum.The bottom panel highlights some gases that were not identified in these data, but that may be present in future observations of other exoplanets.(https://github.com/afeinstein20/wasp39b_niriss_paper/blob/main/scripts/figure3.py)

2
doublet are best explained by [K/O] 0.1-0.5, equivalent to 1-3 solar (see Extended Data Fig.∼ × 8).While the suggested K/O ratio might be a lower limit owing to the photoionization of K at upper atmospheres25 , these observations demonstrate the power of NIRISS-SOSS to obtain simultaneous inferences on novel atmospheric tracers beyond the well-established C/O ratio.

Figure 4 :
Figure 4: Impact of the carbon-to-oxygen ratio (C/O) and metallicity on the JWST-NIRISS spectrum of WASP-39b.Top: Variation of the C/O in the best-fit reference model, while keeping the metallicity, redistribution, and K/O parameters from the reference model the same, and fitting for the cloud parameters and scaled planetary radius to best explain the observations.Under these equilibrium conditions, increasing the C/O results in less H 2 O and more CH 4 , the latter having spectroscopic signatures incompatible with the observations.To mute these incompatible CH 4 features at high C/O, the model requires a higher degree of cloudiness that also mutes any remaining H 2 O features in the spectrum.Bottom: The same exercise as above, but instead we vary the metallicity parameter.The metallicity constraint is driven by the λ > 2μm data; the high-metallicity models ([M/H]> 2) expect larger transit depths than what is seen in the data.The same reference model is plotted as a thick black line in both panels.(https://github.com/afeinstein20/wasp39b_niriss_paper/blob/main/scripts/figure4.py)

3 .
We therefore fix the ∼ ∼ ′′ "soss_transform" parameter to [0, 0, 0], and perform the extraction with a box size of 25 pixels.Any remaining 5 outliers in the resulting spectra are then identified and clipped.Currently, > σ supreme-SPOON does not explicitly treat contamination from zeroth orders of background stars that intersect the trace.Extended Data Figure 2: Comparison of averaged background frames computed for each reduction pipeline.Top: An out-of-transit (OOT) median image from the Stage 2 output files in data numbers per second (DN s -1 ).Middle: Estimated median background frames for four example pipelines: (b) nirHiss, (c) supreme-SPOON, (d) transitspectroscopy, and (e) iraclis.All reduction pipelines use the predefined zodiacal background model provided on the STScI JDox User Documentation website.nirHiss estimates the 0 th order contaminants by taking a smoothed median from the F277W filter integrations.We note that the background frame from supreme-SPOON (c1) for group eight is shown here, and was scaled by a factor of~0.02 to lie on the same scale as the background from the other three pipelines.iraclis subtracts a median per column to remove additional 1/f noise.Bottom: Background subtracted median integrations for each pipeline, plotted in [DN s -1 ] but scaled from -5 to 5 to highlight subtle changes in the background.For these integrations, we define OOT as integrations 1-200 and 400-518.(https://github.com/afeinstein20/wasp39b_niriss_paper/blob/main/scripts/edfigure2.py)

Figure 4 :
White light curves and residuals between the light curve and best-fit exoplanet transit model per each reduction pipeline.Left: White light curves for Order 1 (top) and Order 2 (bottom) with the out-of-transit scatter noted in the figure text.Right: The residuals, in ppm, between the plotted light curves and best-fit exoplanet transit model.The start of transit ingress and end of transit egress are marked with dashed vertical lines; the transit midpoint is marked with a solid vertical line.(https://github.com/afeinstein20/wasp39b_niriss_paper/blob/main/scripts/edfigure4.py) the spectral light curves, fitting only for the , fixing the orbital   / * parameters, and , and the transit mid-time to the above white results.In both cases the / *  models also included a quadratic detrending function that was multiplied by the transit model.After modelling, we applied a spectral decontamination step, taking advantage of the varying total flux across the spectral traces.Due to the contamination we have:, where is the out-of-transit flux (star and contamination) and is  the contaminating source.Hence, for each wavelength we fitted for and applied the  correction: .This procedure is effective in removing  /  −  ( )

2 2 /
specific values per model indicated in the legend of Extended Data Fig.6.All of our forward model grids consistently indicate super-solar metallicity ([M/H]= -) and sub-solar C/O ratio.1Each best-fit spectrum shows different structures at >2 m, as the spectra at these wavelengths µ are more sensitive to the treatment of cloud properties (see next subsection for details).The best-fit spectra from PICASO, ATMO, and Phoenix indicate atmospheric metallicities of [M/H]=1.7,1.0, and 2.0, respectively.These models also consistently indicate C/O ratios between -, corresponding to the lowest C/O ratio grid point in each grid (see the 0. 229 0. 389 main text for why models prefer lower C/O ratios).Thus, the super-solar metallicity and sub-solar C/O ratio of WASP-39b are consistent across the different model interpretations of the NIRISS-SOSS transmission spectrum.We also find that clear atmospheric models fail to fit the observed spectrum even at very high metallicity ([M/H]=2.0),as shown in the bottom panel of Extended Data Fig.6.The clear models fail to match the amplitudes of H 2 O absorption features at , , , and 1.8μm λ = 0. 9 1. 15 1. 4 simultaneously.The clear ATMO models fit the data better than the clear PICASO and Phoenix models because the ATMO grid allows lower heat redistribution factors (i.e., cooler atmosphere).The clear models also overestimate the transit depth at ~2 m because of a strong CO λ µ 2 absorption resulting from the inferred high metallicity ([M/H]=2.0).The inability of clear atmosphere models to fit the overall NIRISS spectrum strongly indicates the presence of clouds in the atmosphere, and emphasises the ability of the NIRISS wavelength coverage to break the cloud property-metallicity degeneracy.The best-fit cloud properties are f sed =1 and K zz = 10 9 cm 2 s -1 for Virga clouds in PICASO, a gray cloud opacity of the H 2 Rayleigh scattering opacity 5 × at 0.35 m for ATMO, and a gray cloud top pressure of 3 10 -4 bar for Phoenix.µ × Extended Data Figure 6: A summary of pre-computed forward model fits to NIRISS-SOSS spectrum.Top: Each coloured line shows the best-fit spectrum from the PICASO, ATMO, and PHOENIX cloudy grid.The chi-squared values per number of data points (N obs = 327) are N obs = 2.98, 3.24, and 3.51 for the PICASO, ATMO, and PHOENIX grids, respectively.All grid models consistently indicate a super-solar metallicity of [M/H] = 1-2 and a sub-solar C/O ratio.Bottom: The same as the top panel, but for the best-fit clear atmosphere models.The clear models yield noticeably worse fits to the data: 2 /N obs = 7.02, 4.11, and 8.55 for the PICASO, ATMO, and PHOENIX grids, respectively, which strongly indicates the presence of clouds in the atmosphere.

×Extended Data Figure 8 :
Using the reference model atmospheric properties, (e.g., [M/H]=1.37,C/O=0.2, full redistribution f=1), we search for the best-fit [K/O] while simultaneously adjusting the 1 bar radius and the parameters for the inhomogeneous cloud model 3, when only using the observations blue-wards of 0.8 m.The best-fit [K/O]=0.4 is consistent with the inferences using µ all the data and the data blue-wards of 2.0 m only.This model is shown in Extended Data µ Figure 8 in green.For the best fit cloud parameters and 1 bar radius, we compute a series of K/O ratios spanning sub-solar and super-solar values.Our results find that sub-solar K/O ratios are disfavoured to , while super-solar values are disfavoured to .Evidence for super-solar [K/O] in WASP-39b.We fit for [K/O] while keeping the rest of the model parameters (e.g., C/O, metallicity, and redistribution) the same as our reference model and fitting for the cloud parameters and scaled planetary radius.Here, we present the different [K/O] models (solid lines) we fit against the transmission spectrum at R=300 (black and white points).We represent each model's respective fit in the orange shading.(https://github.com/afeinstein20/wasp39b_niriss_paper/blob/main/scripts/edfigure8.py) . The line lists of several key species are: H 2 O 76 , CH 4 CH, CH 4 , CN, CO, CO 2 , COF, C 2 , C 2 H 2 , C 2 H 4 , C 2 H 6 , CaH, CrH, FeH, HCN, HCl, HF, HI, HDO, HO 2 , H 2, H 2 S, H 2 O, H 2 O 2, H 3 +, MgH, NH, NH 3 , NO, N 2 , N 2 O, OH, O 2 , O 3 , PH 3 , SF 6 , SiH, SiO, SiO 2 , TiH, TiO, VO, and atoms up to U. The line list of H 2 O is from BT2 76 , other molecular lines from HITRAN 2008 CO, CO 2 , C 2 H 2 , C 2 H 4 , C 2 H 6 , CrH, Cs, Fe, FeH, HCN, H 2 , H 2 O, H 2 S, H 3 +, OCS, K, Li, LiCl, LiH, MgH, NH 3 , N 2 , Na, PH 3 , Rb, SiO, TiO, and VO.The line lists of several key species are: H 2 O 90 , The opacity sources considered are H 2 -H 2 and H 2 -He CIA 107 , H 2 O 90,108 , CO 2

Extended Data Table 1: An outline of reduction and fitting pipelines used to produce transmission spectra for WASP-39b with NIRISS-SOSS.
Size of the box aperture is listed in parentheses when appropriate.All spectra will be made publicly available.

Table 2 : White-light curve best-fit orbital parameters from Order 1.
Transit time (t 0 ) is presented with respect to t 0 -2459787 [BJD].