JWST transmission spectroscopy of HD 209458b: a super-solar metallicity, a very low C/O, and no evidence of CH₄, HCN, or C₂H₂

Detected initially by the radial velocity method in 1999, the hot Jupiter HD 209458b was the first exoplanet found to transit its host star (Charbonneau et al., 2000; Henry et al., 2000). Since then, it has been one of the most frequently studied exoplanets and it has been the subject of a number of breakthroughs that sparked the study of exoplanetary atmospheres. Via transmission spectroscopy, it was the subject of the first exoplanet atmosphere detection (Charbonneau et al., 2002), the first found to have an escaping atmosphere (Vidal-Madjar et al., 2003), the first observed to possess atomic carbon and oxygen (Vidal-Madjar et al., 2004), and the first to have its orbital velocity measured, thus turning the system into a double-lined spectroscopic binary (Snellen et al., 2010). It was one of the first two exoplanets with detected infrared emission (Deming et al., 2005), and it was one of the first two exoplanets with the spectroscopic identification of water (Deming et al., 2013a). What’s more, it was once thought to possess a stratospheric temperature inversion (Knutson et al., 2008; Burrows et al., 2007), although that hypothesis was later refuted by subsequent observations and analysis (Diamond-Lowe et al., 2014; Schwarz et al., 2015; Line et al., 2016).

HD 209458b is still among the best targets for atmospheric study in this era of thousands of known transiting planets. It has a relatively bright host star, a favorable planet-to-star radius ratio, a high equilibrium temperature, and a low surface gravity. Ultimately, it has the highest transmission spectroscopy metric ($\sim 900$¹¹1https://tess.mit.edu/science/tess-acwg/) of all known exoplanets (Kempton et al., 2018).

Despite its great observability, several fundamental questions about the composition of HD 209458b’s atmosphere remain unsolved. One of these questions is its atmospheric water abundance. Several groups have derived sub-solar water abundances from measurements of its atmosphere (Madhusudhan et al., 2014; MacDonald & Madhusudhan, 2017a; Brogi et al., 2017; Welbanks et al., 2019a; Pinhas et al., 2019), with a recent re-analysis providing estimates consistent with solar expectations within the uncertainties (Welbanks & Madhusudhan, 2021). This could be caused by an overall low metallicity or just a low oxygen abundance in the planet, neither of which are expected by standard models for giant planet formation (Öberg & Bergin, 2016; Booth et al., 2017). However, Line et al. (2016), Sing et al. (2016) and Tsiaras et al. (2018) also reported solar to super-solar H₂O abundances, which is more in line with expectations from traditional planet formation models (Owen & Encrenaz, 2006; Öberg et al., 2011).

A second open question surrounds the carbon and nitrogen chemistry. Interpretation of the Hubble Space Telescope (HST) transmission spectrum initially suggested strong evidence for NH₃ and/or HCN (MacDonald & Madhusudhan, 2017a), but a subsequent analysis lowered the detection significance and highlighted NH₃ as being the more likely of the two (MacDonald & Madhusudhan, 2017b). One the other hand, two high-resolution spectroscopy (HRS) analyses both claimed the presence of HCN (Hawker et al., 2018; Giacobbe et al., 2021). The latter of these two studies also claimed the detection of NH₃, CH₄, and C₂H₂ in addition to H₂O and CO (Giacobbe et al., 2021). When assuming equilibrium chemistry and a clear atmosphere, the presence of all of these molecules together suggests a highly sub-solar metallicity ($<$1% of solar) and/or a C/O ratio of around 1 or higher (compare to the solar value of 0.59, Asplund et al., 2021), which together challenge planet formation models (Mousis et al., 2012; Madhusudhan et al., 2014).

In this paper, we present the first transmission spectrum of the archetypical hot Jupiter HD 209458b obtained with the James Webb Space Telescope (JWST) to answer these longstanding questions. We describe our observations and data analysis in §2, atmospheric modeling in §3 and results in §4.

We observed two transits of HD 209458b with JWST’s NIRCam instrument (Greene et al., 2017) on November 10 and 16, 2022 (program GTO 1274, J. Lunine PI). Each observation lasted 8.01 hours, which is long enough to sample the 3.12-hour transit and the baseline flux before and after the transit. Both observations used the module A grism R mode to obtain time-series near-infrared spectra in the long-wavelength (LW) channel. The first observation collected data with the F444W filter, yielding spectra from 3.86 to 5.06 $\mu m$. The second observation used the F322W2 filter, yielding spectra from 2.36 to 4.02 $\mu m$.

Both visits employed the SUBGRISM64 subarray and BRIGHT2 readout pattern. The first observation used six groups per integration for a total of 7,107 integrations, while the second used four groups per integration for a total of 9,473 integrations. Photometry was also obtained simultaneously in the short-wavelength (SW) channel during each LW channel observation. The SW photometry for both visits was obtained in a narrow band at 2.12 $\mu m$ using the WLP4 filter with the BRIGHT1 readout pattern.

We reduced and analyzed the LW data independently using two different pipelines. The first pipeline we used is Eureka! (Bell et al., 2022), which has been utilized extensively by the JWST Transiting Exoplanet Early Release Science (JTEC ERS) Team (JTEC ERS Team et al., 2023; Ahrer et al., 2023; Alderson et al., 2023; Feinstein et al., 2023; Rustamkulov et al., 2023; Coulombe et al., 2023). The second pipeline we used is SPARTA²²2https://github.com/ideasrule/sparta, which was first utilized for the MIRI/LRS phase curve of GJ 1214b (Kempton et al., 2023) and is also now being included in the JTEC ERS analyses of MIRI/LRS data of WASP-43b (Bell et al., submitted) and WASP-39b (Powell et al., submitted).

Our Eureka! implementation follows the analyses of NIRCam data presented in Ahrer et al. (2023) and Bean et al. (2023). The first two stages in Eureka! are identical with those in the JWST Science Calibration Pipeline (jwst) except that we did a group-level background subtraction and we increased the jump detection threshold to 6.0. Following Stage 2, Eureka! performed a column-by-column linear fit to calculate the background in the region beyond 13 pixels relative to the center of the spectral trace and subtracted it from the region of interest. Correction of the curvature of the spectral trace was performed by shifting columns to align the center of the spectral trace along the same row. Then optimal spectral extraction as defined in Horne (1986) was performed for each integration within eight pixels on both sides of the trace. The spatial profile used in the optimal extraction weighting is based on a median frame cleared of 10$\sigma$ outliers. We excluded the five spectroscopic light curves with scattering factors greater than 1.6$\times$ photon noise in the presented spectra.

SPARTA is a new, end-to-end pipeline that begins with the raw, uncalibrated files with JWST. SPARTA has its own up-the-ramp fitting that does not depend on the JWST pipeline. The detailed reduction algorithms of SPARTA can be found in Kempton et al. (2023). For F332W2 and F444W, we let the spectral center be at the 35th and 33rd pixel respectively and the extraction window to be 8 pixels. For spectroscopic light curve fitting, we excluded three data points with scattering factors greater than 1.35 or smaller than 1.0. As with Eureka!, the transit light curve parameters were estimated using the dynamic nested sampling algorithm (Higson et al., 2018) as implemented by the dynesty package³³3https://dynesty.readthedocs.io/en/stable/.

We generated both “white” light curves that were summed over the full bandpass of each observation and spectroscopic light curves that were summed over 0.04 $\mu m$ each (yielding 69 channels in total) for the LW data. A 9$\sigma$ outlier rejection on each light curve were performed and we fit each light curve with a transit model from the batman code (Kreidberg, 2015) combined with a systematics model that is linear with time (i.e., $c_{0}+c_{1}t$). We trimmed the first 990 integrations of the first observation and 1,176 integrations of the second (both approximately the first 1 hour) due to the strong exponential-like ramp at the beginning of the observations.

We adopted the orbital period as 3.52474859 days (Stassun et al., 2017), eccentricity as 0.0, and argument of periapsis as 90.0$\arcdeg$. The inclination $i$ and semi-major axis in units of the host star radius $a/R_{s}$ were determined by fitting the white light curve of the second visit because it has less noise than the first visit (more photons were collected in these shorter wavelength data). Then $i$ and $a/R_{s}$ were fixed in the analysis of all the spectroscopic channels. The mid-transit time of the two observations is determined by fitting the white light curve of each observation. The best estimated parameters can be found in Table 1.

As is typical, the measured transit depths depend on the limb darkening assumptions. We tried using limb-darkening coefficients for a four-parameter “non-linear” law calculated from 3D stellar model atmospheres specific to HD 209458 (Hayek et al., 2012), but the results do not match our data well. Therefore, we elected to instead fit for the limb darkening coefficients in the light curve modeling. We adopted the quadratic limb-darkening law re-parameterised by Kipping (2013).

The SW data were reduced by enabling photometric analysis in Eureka! (see Bean et al., 2023). First, the centroid of the image was determined by a 2D Gaussian fit. We then performed aperture photometry with radius of 45 pixels and a background annulus spanning from 100 to 120 pixels because this combination minimized the scatter in the light curve. We didn’t correct for $1/f$ noise because it was not evident in our data. As we did for the LW data, the light curves from SW were fitted using dynamic nested sampling.

The transmission spectrum of HD 209458b measured in the LW and SW NIRCam data by Eureka! is shown in Figure 1. The weighted mean difference between the independent Eureka! and SPARTA reductions of the LW data is 1.67 $\sigma$ (see Figure 6 for the comparison). The overlapping region of the two LW filters (3.90 to 4.00 $\mu m$) agrees at 1.9 $\sigma$ for Eureka! reduction and 0.37$\sigma$ for SPARTA. However, we didn’t include the SW data in the atmospheric retrieval (see next section) because of the high scatter seen in the light curves ($\sim$ 10$\times$ photon noise).

We retrieved the atmospheric properties of HD 209458b by fitting the JWST spectrum using PLATON⁴⁴4https://platon.readthedocs.io/en/latest/ (Zhang et al., 2019, 2020), which has been used to determine the properties of several hot Jupiters (Jiang et al., 2021; Ahrer et al., 2022; Spyratos et al., 2023; Bean et al., 2023; August et al., 2023). Our retrieval setup assumed an isothermal atmosphere with equilibrium chemistry and a “patchy” grey cloud deck (the latter motivated by Line & Parmentier, 2016). The retrieval included the planet radius at 10⁵ Pa, temperature, metallicity ([M/H] = log(M/H) - log(M/H)_Sun), carbon-to-oxygen (C/O) ratio, cloud-top pressure, and cloud fraction on the day-nightside terminator as free parameters. For the cloud fraction we use the prescription of Pinhas et al. (2019). Additionally, we included the mixing ratios of CH₄, NH₃, C₂H₂, and HCN as free parameters (with log-uniform priors from $10^{-10}$ to $10^{-2}$) to obtain constraints on their abundances separate from the equilibrium chemistry predictions because these four species are key detections in Giacobbe et al. (2021).

Two additional cloud models are tested as follows. The first one is with spectral slope (where the absorption coefficient is characterized by ${\alpha\propto A\times\lambda^{slope}}$)(Zhang et al., 2019). Besides the six parameters described above, we retrieved scattering amplitude $A$ and slope, but neither of them can be well constrained. The second is Mie scattering which has a complex refractive index of ${1.33-0.1j}$. We retrieved number density and the size of particles and the resulting C/O and [M/H] are consistent with the simpler “patchy” grey cloud deck model described above.

The abundance grid is computed by FastChem⁵⁵5https://github.com/exoclime/FastChem and has five dimensions: species name (in total 34 atomic and molecular species, same as Zhang et al. (2019)), temperature (100 – 3000K with 100K interval), pressure($10^{-4}$ – $10^{8}$ Pa with decade interval), metallicity(from [M/H] = -1 to [M/H] = 2 with 0.03 interval) and C/O (0.001 to 2.0⁶⁶6the interval is non-linear: [0.001, 0.005, 0.01, 0.02, 0.03, 0.04, 0.05, 0.1, 0.2, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 1.0, 1.05, 1.1, 1.2, 1.4, 1.6, 1.8, 2.0]). The carbon abundance is computed by scaling the solar oxygen abundance with different C/O, and the other elements’ abundances (except H and He) are scaled with metallicity.

The absorption coefficients used in the modeling are from the DACE opacity database⁷⁷7https://dace.unige.ch. with a resolution of $R=100,000$ over 0.5 – 12 $\mu$m. We included CH₄, CO, CO₂, H₂O, H₂S, NH₃, C₂H₂, HCN, and SO₂ as these are the molecules making major contribution within this wavelength range (Pinhas et al., 2019). The adopted absorption coefficients that are compatible with PLATON used in our retrieval can be found in §6.

As a consistency check, we used PLATON to retrieve the properties of WASP-39b from the JTEC ERS NIRSpec G395H spectrum of the planet (Alderson et al., 2023). We removed the data from 3.9 $\mu$m to 4.1 $\mu$m on WASP-39b’s spectrum that covers the SO₂ absorption feature produced by photochemistry (Tsai et al., 2023). The best fit has T = 982${}^{+40}_{-39}$ K, [M/H] = 1.39$\pm 0.16$ and C/O = 0.66${}^{+0.11}_{-0.25}$, which are consistent with the results in Alderson et al. (2023) and Constantinou et al. (2023). However, we found a strong degeneracy between temperature and $R_{p}$ (see more discussion of this in the next section).

In Figure 1(a) we show the best-fit model to our JWST transmission spectrum of HD 209458b, with the contributions from H₂O, CO₂, CO, H₂S and the patchy cloud highlighted. The absorption features in the spectrum are primarily due to water (feature centered at 2.8 $\mu m$) and carbon dioxide (feature centered at 4.3 $\mu m$), with perhaps a minor contribution from H₂S. The water and carbon dioxide features are reduced in amplitude by about a factor of two due to the presence of a cloud deck. Our data favor a cloud patchiness fraction of $\sim$68% at 3.0$\sigma$ significance⁸⁸8calculated as median divided by deviation. Previous HST/WFC3 observations taken at shorter wavelengths (Deming et al., 2013b) identified water vapour. This is the first detection of CO₂ in HD209458b’s atmosphere, continuing the trend of detections of this molecule using JWST after WASP-39b (JTEC ERS Team et al., 2023), HD 149026b (Bean et al., 2023), and K2-18b (Madhusudhan et al., 2023). There is no evidence for additional absorbers.

In Figure 1(b) we present models where we scale our equilibrium abundance of CH₄, NH₃, C₂H₂, and HCN to the notional abundances from Giacobbe et al. (2021, their Extendend Data Table 4), which gave the maximum cross-correlation signal in their data on a species-by-species basis. These molecules have absorption features in our bandpass and would have shown up (to varying degree) if their abundances were as high as those suggested by Giacobbe et al. (2021). By including the volume mixing ratios of CH₄, NH₃, C₂H₂, and HCN as free parameters in our retrieval, we provide $3\sigma$ upper limits of (see Figure 2) log($\chi_{\mathrm{CH}_{4}}$) = -5.6,  log($\chi_{\mathrm{NH}_{3}}$) = -4.2, log($\chi_{\mathrm{C}_{2}\mathrm{H}_{2}}$) = -5.7, and log($\chi_{\mathrm{HCN}}$) = -5.1. The posteriors for the abundances of all four molecules are consistent with the chemical equilibrium prediction from our best-fit model.

Of the four extra molecules that we explored, only the notional abundance of NH₃ from Giacobbe et al. (2021) is potentially consistent with our retrieval results i.e., our 3$\sigma$ upper limit is within an order of magnitude of their abundance). Our posterior for NH₃ is not bounded on the low end, which is consistent with the constraints for this molecule from both MacDonald & Madhusudhan (2017a) and MacDonald & Madhusudhan (2017b), which are based on HST/WFC3 data. Therefore, the current space-based data are unable to determine if the molecule is present at equilibrium or higher abundances.

To test if the cloud prescriptions would influence the detection of CH₄, NH₃, C₂H₂, and HCN, we repeated the two cloud models described in §3 (paragraph 2) and recomputed the abundances of these four key molecules. We report the constrained $3\sigma$ upper limits of CH₄, C₂H₂, and HCN are all less than $10^{-5}$.

The abundances suggested by Giacobbe et al. (2021) for the other three molecules we explored (CH₄, HCN, and C₂H₂) are highly inconsistent with our results (i.e., our 3$\sigma$ upper limits are at least an order of magnitude lower than their abundances). However, the reported volume mixing ratios from Giacobbe et al. (2021) are from models that maximize the cross-correlation functions for each species individually, and thus are not retrieved values with proper uncertainties. Therefore, the spectra shown in Figure 1(b) might not represent the actual inference from their data, and further analysis is needed to assess the level of agreement.

In Figure 3, we show the corner plot for our chemical equilibrium retrieval on the Eureka! reduction. Similar to our retrieval on WASP-39b, we found a strong correlation between the temperature and the planet radius, which we believe is caused by the limitation of wavelength coverage. Specifically, the spectrum lacks the continuum fully outside molecular absorption bands and the multiple bands of the same molecule that are helpful for breaking degeneracies in transmission spectra (Benneke & Seager, 2012). On the other hand, these JWST data precisely resolve the shape of the H₂O and CO₂ features, and our assumption of chemical equilibrium provides additional constraints on the retrieval.

To test the robustness of the wavelength coverage of our JWST data, we compared it with the best-fit parameters from joint HST/WFC3 and JWST fitting, and the retrieved [M/H], C/O are consistent at greater than 0.5$\sigma$. We also conducted a test to see the impact of the minor difference in the spectrum obtained by a different data reduction pipeline. We found the retrieved parameters from Eureka!’s spectrum are within 1$\sigma$ of those from SPARTA’s spectrum. The best retrieved atmospheric properties can be found in Table 2.

Using only the JWST data, PLATON favors a metallicity of $3^{+4}_{-1}$ $\times$ solar and a C/O = $0.11^{+0.12}_{-0.06}$ (bolded column in Table 2). The observed planetary atmospheric metallicity-mass trend in our solar system has motivated a number of studies (Thorngren et al., 2016; Welbanks et al., 2019b). As can be seen in Figure 4, our retrieved M/H is within $1\sigma$ of the trend of methane abundances in the solar system giant planets. Our findings suggest that HD 209458b might have undergone a similar amount of planetesimal accretion as the solar system giant planets.

The interior model in Thorngren & Fortney (2019) anticipates the maximum metal enrichment for an exoplanet population given its mass and radius for a “core-less” planet (i.e., the metals and gas are thoroughly mixed within the whole planet). These upper limits (as shown with blue circles in Figure 4) greatly exceed the observed Jupiter and Saturn metal enrichments, implying that a huge percentage of metals must be trapped within a core. The atmospheric metallicity of HD 209458b below this limit indicates that some of the accreted metal is also in its core.

The constraint on the atmospheric C/O of HD 209458b from our spectrum is driven by a combination of the detected H₂O and CO₂, and the lack of detection of other molecules like CH₄ that would be present in atmospheres with higher C/O values. In Figure 1(b), we show models with C/O = 0.5, 1.0, and 1.5. We find that with a higher C/O ratio, the water feature at 2.3 $\mu m$ is weakened and the CO₂ feature is enhanced until both of them are absent for C/O $>$ 1. CH₄ takes up the high carbon abundance above C/O ratios of unity.

The low C/O ratio for HD 209458b inferred in our work is due to the strong H₂O absorption feature, which indicates a high oxygen abundance. As a comparison, we plotted the JWST transmission spectra in units of scale height for WASP-39b and HD 209458b in Figure 5 (left). HD 209458b exhibits a relatively stronger H₂O feature and a weaker CO₂ feature in comparison to WASP-39b. The latter exoplanet has a C/O ranging from 0.3 to 0.46, as reported in the NIRSpec G395H paper (Alderson et al., 2023). The C/O from this particular paper is chosen due to its bandpass similarity to our work. In Figure 5 (right), we show the ratio of H₂O and CO₂ abundances as a function of C/O ratio expected from chemical equilibrium. The $\chi_{\mathrm{H}_{2}\mathrm{O}}/\chi_{\mathrm{CO}_{2}}$ has a dependence on metallicity because CO₂ itself is a strong function of metallicity. Nevertheless, given the similar metallicities of the two planets (WASP-39b is $\sim$ 10 $\times$ solar, JTEC ERS Team et al., 2023), the $\chi_{\mathrm{H}_{2}\mathrm{O}}/\chi_{\mathrm{CO}_{2}}$ ratio is mostly indicative of the different C/O ratios of the planets. The relative strength of H₂O vs. CO₂ absorption thus demonstrates that HD 209458b’s C/O ratio is significantly lower than that of WASP-39b.

In addition to the main atmospheric retrieval described above, we performed a grid fit using the ScCHIMERA Radiative Convective Equilibrium (RCE) solver first described in Piskorz et al. (2018) and recently used in JWST observations of WASP-39b (Rustamkulov et al., 2023), WASP-96b (Radica et al., 2023), and WASP-80b (Bell et al., 2023). We computed a grid of models under the assumption of 1-dimensional RCE given an irradiation and elemental composition. The grid is computed for steps in the irradiation temperature of $T_{\rm irr}$ (1347 – 1557 K in steps of 15 K), [M/H] (-0.5 – 2.25 in steps of 0.125), and C/O (0.1 – 0.75 in steps of 0.05). A detailed description of the ScCHIMERA grid and parameter estimation is available in Radica et al. (2023) and Bell et al. (2023).

Generally, we performed the parameter estimation over the ScCHIMERA grid by post-processing the 1D-RCE atmospheric structures through a transmission spectrum routine while considering the presence of inhomogeneous clouds and power-law hazes. The resulting parameter estimation derived a metallicity of $\sim 1-2\times$ solar and a sub-solar C/O ratio. The C/O ratio runs up against the lower bound limit of the grid (0.1) and has a 2$\sigma$ upper limit of 0.32. Given the consistency of the results from the two different types of retrievals we emphasize the PLATON results as the main finding in this paper.

Kawashima & Min (2021) report large differences in CH₄ abundance and C/O between equilibrium and disequilibrium limits when retrieving spectra covering 2.5 – 4.0 $\mu$m. By introducing eddy diffusion transport, they found CH₄ quenches at P = 1 bar thus resulting in less CH₄ compared to the equilibrium case and C/O is raised to 0.5 in disequilibrium as compared to 0.19 in equilibrium. In order to test the impact of chemical disequilibrium on C/O, on top of the six parameters discussed in §3, we added one more free parameter P_quenching in PLATON, the pressure at which quenching happens. However, P_quenching cannot be constrained according to our retrieval. Though our results disfavor the disequilibrium scenario, this may partly be because we assume the same quenching pressure for all molecules, while in reality it should be different (Moses, 2014). We also assume an isothermal atmosphere with fixed abundance profile while in Kawashima & Min (2021) the profile varies within the retrieval. Further modeling work on this data set using more sophisticated treatments of disequilibrium chemistry are thus warranted.

In this paper, we present the transmission spectrum of the transiting hot Jupiter HD 209458b with data observed with JWST/NIRCam from 2.3 – 5.1 $\mu$m. The data show clear features from water, carbon dioxide, and clouds. We do not detect evidence for additional molecules that had been previously claimed for this planet. Our retrieval results suggest a mild atmospheric metallicity enhancement between that of Jupiter and Saturn in our own solar system. The data also suggest a very low C/O ratio that stands out from other emerging JWST results for giant exoplanets.

In terms of non-detections, our upper limits on the abundances of CH₄, HCN, and C₂H₂ in particular make the presence of these molecules in the atmosphere of HD 209458b as claimed by Hawker et al. (2018, HCN only) and Giacobbe et al. (2021, all three molecules) controversial. Nevertheless, it is important to point out that Giacobbe et al. (2021) did not perform a retrieval to put formal constraints on the abundances of the molecules they detected due to the challenge of such analyses on HRS data (e.g., Brogi & Line, 2019). Therefore, it is not clear what the statistical significance is of the possible tension between the results. HRS in principle may be sensitive to trace species that elude low-resolution spectroscopy, but it remains to be seen whether the very low abundances constrained by our data would still yield a detection in HRS data.

Studies on the robustness of molecular detections by HRS of exoplanets demonstrate that some detrending methods may induce false positive or inflated detections (Cheverall et al., 2023). In their case study on HD 209458b, HCN, NH₃ and CH₄ were not observed by Cheverall et al. (2023). However, a signal for CH₄ can be detected if the telluric contamination is not correctly removed. While methane (CH₄) is ubiquitous in the atmosphere of solar system giant planets, it absence has long been one of the core puzzles of the study of exoplanetary atmospheres (Gibson et al., 2011; Benneke et al., 2019; Baxter et al., 2021; JTEC ERS Team et al., 2023). However, recent JWST observations on transiting exoplanets WASP-80 b and K2-18 b show evidence of methane (Bell et al., 2023; Madhusudhan et al., 2023). Understanding this molecule’s absence or existence is essential for understanding how planets form and evolve, how atmospheric processes work, and the habitability of planets. Future JWST observations, like cycle 2 program 3557, might aid in solving this mystery.

Our study disproves the previous claims of low water abundance for this planet (Madhusudhan et al., 2014; MacDonald & Madhusudhan, 2017a; Brogi et al., 2017; Welbanks et al., 2019a; Pinhas et al., 2019). Instead of being depleted in water, our best-fit model gives enhanced metallicity and low C/O, indicating the atmosphere of HD 209458b is rich in oxygen. We found a metallicity that is consistent with, but more precise than Line et al. (2016) in emission and Welbanks & Madhusudhan (2021) in transmission, which extends the agreement that is seen between the solar system trend and exoplanet atmosphere abundances.

Line et al. (2016) and Brogi et al. (2017) have constrained the C/O of HD 209458b to $<1$, yet our very low C/O ratio provides valuable insights into the planetary formation and evolution. Espinoza et al. (2017) suggests a low C/O in gas giants compared to parent stars is caused by metal enrichment but not dependent on the formation location. Additionally, planetesimal pollution (Öberg et al., 2011) caused by formation and migration sufficiently inward of the snowlines of carbon-bearing species may also result in low ($<0.5$) C/O and elevated metal enrichment. Through a comparative analysis of the spectra of this study and WASP-39b, we found the two have similar molecular features (namely water and CO₂), while the variation in the ratio of water to carbon dioxide abundance leads to a distinct difference in the C/O. This is primarily because the value of $\chi_{\mathrm{H}_{2}\mathrm{O}}/\chi_{\mathrm{CO}_{2}}$ serves as a robust indicator of C/O, assuming comparable metallicity. It is important to note that this relationship is not influenced by specific retrieval models, thus revealing the intrinsic characteristics of such planets.

Tsai et al. (2023) showed evidence of photochemically-produced SO₂ in the atmosphere of WASP-39b. Since we do not detect SO₂ at 4.05 $\mu$m, we have not made an effort to study the potential impact of chemical disequilibrium on the H₂S and SO₂ abundance. Detailed modeling to assesss our JWST transmission spectrum in the context of the disequilibrium process would be interesting.

The data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute. The specific observations analyzed can be accessed via https://doi.org/10.17909/f5j3-jq48 (catalog DOI: 10.17909/f5j3-jq48). The data that were used to create all of the figures will be freely available on Zenodo (Xue et al., 2024). All additional data is available upon request.