With sizable volatile envelopes but smaller radii than the solar system ice giants subNeptunes have been revealed as one of the most common types of planet in the galaxy. While the spectroscopic characterization of larger subNeptunes 2.54 R SUBSUB has revealed hydrogendominated atmospheres smaller subNeptunes 1.62.5 R SUBSUB could either host thin rapidly evaporating hydrogenrich atmospheres or be stable metalrich water worlds with high mean molecular weight atmospheres and a fundamentally different formation and evolutionary history. Here we present the 0.62.8 m JWSTNIRISSSOSS transmission spectrum of GJ 9827 d the smallest 1.98 R SUBSUB warm T SUBeqA0.3SUB  620 K subNeptune where atmospheric absorbers have been detected to date. Our two transit observations with NIRISSSOSS combined with the existing HSTWFC3 spectrum enable us to break the cloudsmetallicity degeneracy. We detect water in a highly metalenriched steam world atmosphere OH of 4 by mass and HSUB2SUBO found to be the background gas with a volume mixing ratio of gt31. We further show that these results are robust to stellar contamination through the transit light source effect. We do not detect escaping metastable He which combined with previous nondetections of escaping He and H supports the steam atmosphere scenario. In waterrich atmospheres hydrogen loss driven by water photolysis happens predominantly in the ionized form which eludes observational constraints. We also detect several flares in the NIRISSSOSS light curves with farUV energies of the order of 10SUP30SUP erg highlighting the active nature of the star. Further atmospheric characterization of GJ 9827 d probing carbon or sulfur species could reveal the origin of its high metal enrichment.
['Exoplanet atmospheres', 'Exoplanet atmospheric composition', 'Exoplanet atmospheric evolution', 'Exoplanet structure', 'Planetary atmospheres', 'Exoplanet astronomy', '487', '2021', '2308', '495', '1244', '486', 'Astrophysics - Earth and Planetary Astrophysics', 'Astrophysics - Solar and Stellar Astrophysics']
JWSTNIRISS Reveals the Waterrich Steam World Atmosphere of GJ 9827 d
{"JWST/NIRISS reveals the water-rich 'steam world' atmosphere of GJ 9827 d": "Caroline Piaulet-Ghorayeb , 1 Bjorn Benneke , 1 Michael Radica , 1 Eshan Raul , 2 Louis-Philippe Coulombe , 1 Eva-Maria Ahrer , 3 Daria Kubyshkina , 4 Ward S. Howard , 5 Joshua Krissansen-Totton , 6 Ryan MacDonald , 2 Pierre-Alexis Roy , 1 Amy Louca , 7, 8 Duncan Christie , 3 Marylou Fournier-Tondreau , 9 Romain Allart , 1 Yamila Miguel , 7, 8 Hilke E. Schlichting , 10 Luis Welbanks , 11 Charles Cadieux , 1 Caroline Dorn , 12 Thomas M. Evans-Soma , 13, 3 Jonathan J. Fortney , 14 Raymond Pierrehumbert , 9 David Lafreni'ere , 1 Lorena Acu˜na , 3 Thaddeus Komacek , 15 Hamish Innes , 16, 17 Thomas G. Beatty , 18 Ryan Cloutier , 19 Ren'e Doyon , 1 Anna Gagnebin , 14 Cyril Gapp , 3 and Heather A. Knutson 20  \n1 Department of Physics and Trottier Institute for Research on Exoplanets, Universit'e de Montr'eal, Montreal, QC, Canada caroline.piaulet@umontreal.ca  \n2 Department of Astronomy, University of Michigan, Ann Arbor, MI, USA 3 Max Planck Institute for Astronomy, Konigstuhl 17, D-69117 Heidelberg, Germany 4 Space Research Institute, Austrian Academy of Sciences, Schmiedlstrasse 6, 8042 Graz, Austria 5 Department of Astrophysical and Planetary Sciences, University of Colorado, 2000 Colorado Avenue, Boulder, CO 80309, USA 6 Department of Earth and Space Sciences/Astrobiology Program, University of Washington, Seattle, WA, USA 7 Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands 8 SRON Netherlands Institute for Space Research, Niels Bohrweg 4, 2333 CA Leiden, The Netherlands 9 University of Oxford, Department of Physics Oxford, OX1 3PW, UK 10 Department of Earth, Planetary, and Space Sciences, University of California, Los Angeles, Los Angeles, CA 90095, USA 11 School of Earth & Space Exploration, Arizona State University, Tempe, AZ 85257, USA 12 Department of Physics and Institute for Particle Physics and Astrophysics, ETH Zurich, CH-8093 Zurich, Switzerland 13 School of Information and Physical Sciences, University of Newcastle, Callaghan, NSW, Australia 14 Department of Astronomy & Astrophysics, University of California, Santa Cruz, CA 95064, USA 15 Department of Astronomy, University of Maryland, College Park, MD 20742, USA 16 Department of Earth Sciences, Freie Universitat Berlin, Malteserstr. 74-100, 12249 Berlin, Germany 17 Institute of Planetary Research, German Aerospace Center (DLR), Rutherfordstraße 2, 12489 Berlin, Germany 18 Department of Astronomy, University of Wisconsin-Madison, Madison, WI 53706, USA 19 Department of Physics & Astronomy, McMaster University, 1280 Main Street W, Hamilton, ON L8S 4L8, Canada 20 Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, CA 91125, USA", 'ABSTRACT': "With sizable volatile envelopes but smaller radii than the solar system ice giants, sub-Neptunes have been revealed as one of the most common types of planet in the galaxy. While the spectroscopic characterization of larger sub-Neptunes (2.5-4R ⊕ ) has revealed hydrogen-dominated atmospheres, smaller sub-Neptunes (1.6-2.5R ⊕ ) could either host thin, rapidly evaporating hydrogen-rich atmospheres or be stable metal-rich 'water worlds' with high mean molecular weight atmospheres and a fundamentally different formation and evolutionary history. Here, we present the 0.6-2.8 µ m JWST NIRISS/SOSS transmission spectrum of GJ 9827 d, the smallest (1.98 R ⊕ ) warm (T eq , A B =0 . 3 ∼ 620 K) sub-Neptune where atmospheric absorbers have been detected to date. Our two transit observations with NIRISS/SOSS, combined with the existing HST /WFC3 spectrum, enable us to break the cloudsmetallicity degeneracy. We detect water in a highly metal-enriched 'steam world' atmosphere (O/H of ∼ 4 by mass and H 2 O found to be the background gas with a volume mixing ratio (VMR) of > 31%). We further show that these results are robust to stellar contamination through the transit light source effect. We do not detect escaping metastable He, which, combined with previous nondetections of escaping He and H, supports the steam atmosphere scenario. In water-rich atmospheres, hydrogen loss driven by water photolysis happens predominantly in the ionized form which eludes observational constraints. We also detect several flares in the NIRISS/SOSS light-curves with far-UV energies of the order of 10 30 erg, highlighting the active nature of the star. Further atmospheric characterization of GJ 9827 d probing carbon or sulfur species could reveal the origin of its high metal enrichment.  \nKeywords: Exoplanets (498); Exoplanet atmospheres (487); Planetary atmospheres (1244)", '1.1. The radius valley': "One of the most fundamental observational results from the study of exoplanets has been the discovery that small, close-in planets are bifurcated into two seemingly separate populations (Fulton et al. 2017). The smallplanet size distribution displays a feature known as the 'radius valley': an observed dearth of close-in small planets with sizes of 1.8 to 2.0 Earth radii around FGK stars. This statistical feature has been proposed to be a result of photoevaporation (Owen & Wu 2017) where the mass loss is driven by the high-energy irradiation from the host star, core-powered mass loss (Ginzburg et al. 2018; Gupta & Schlichting 2020) whereby hydrogen envelopes are lost as a result of the bolometric luminosity from the host star and the residual internal heat which is released from the planetary core as it cools down over time; or even the superposition of gas-rich and gas-poor planets that form before and after disk dissipation (Lee & Connors 2021; Lopez & Rice 2018; Zeng et al. 2019). Recent studies suggest that for lower-mass stars, the gap in the radius distribution becomes more and more populated, which suggests enhanced diversity in the bulk volatile content of small planets orbiting low-mass stars (Luque & Pall'e 2022; Ho et al. 2024; Venturini et al. 2024).", '1.2. What masses tell us about compositions': "In each of these theoretical frameworks, the abundant 'sub-Neptunes' comprising the ≳ 2 R ⊕ population are low-density planets with cores made of a mix of ices and rocky material underlying a hydrogen-dominated envelope. Meanwhile, the smaller 'super-Earths' are understood as rocky in composition, either due to their formation or as a result of atmosphere stripping. The precisely measured bulk densities of the smallest known ( ≲ 1 . 6 R ⊕ ) super-Earths agree well with rock/iron core compositions (Parc et al. 2024), in line with theoretical expectations. The low densities of sub-Neptunes, on the other hand, can be equally well explained by a wider range of compositions with different ratios of H 2 compared to high mean molecular weight (HMMW) volatiles. Plausible compositions compatible with subNeptune densities span the continuum from the canonical view of percent-by-mass H/He envelopes of 'gas dwarfs' (Fortney et al. 2013; Benneke et al. 2024), to 'metal-rich' or 'steam world' atmospheres rich in HMMW volatiles amounting to tens of percent of the  \nplanet mass (Figure 1, see also Aguichine et al. 2021, Piaulet et al. 2023).  \nFurthermore, recent theoretical (Madhusudhan et al. 2021; Innes et al. 2023; Dorn & Lichtenberg 2021) and observational findings (Madhusudhan et al. 2023; Holmberg & Madhusudhan 2024; Benneke et al. 2024; Luque & Pall'e 2022) support greater compositional diversity within the sub-Neptune population. While the coolest planets may be 'Hycean worlds' with water in condensed liquid/ice phases under a H/He atmosphere, the warmer temperature of the vast majority of planets would instead imply a supercritical state for the water/volatile layer. Contrary to liquid water, supercritical water is highly miscible with hydrogen. Such warm, volatile-rich planets may therefore host exposed 'mixed' envelopes where both H 2 and potentially large amounts of HMMW volatiles coexist (Wogan et al. 2024; Burn et al. 2024; Benneke et al. 2024). In addition, water is expected to partition between the planetary envelope, the molten or solid mantle, and the metal core, leading to more uncertainties in planetary structure models (Schlichting & Young 2021; Luo et al. 2024; Dorn & Lichtenberg 2021).  \nMass and radius measurements for individual planets are insufficient to resolve this compositional ambiguity due to an inherent degeneracy in the interior models (Rogers & Seager 2010) which prevents the disentangling of thin unstable hydrogen layers from deep layers of HMMW volatiles. Moreover, this question cannot yet be resolved at the population level (Rogers et al. 2023) when accounting for the effects of both thermal evolution and photoevaporation.", '1.3. Exploring the diverse population of small sub-Neptunes': "Spectroscopic observations of planetary atmospheres can resolve these mass-radius degeneracies by providing precise measurements of the mean molecular weight in the upper atmosphere, modulo our understanding of how the amount of volatiles in the atmosphere translates to a bulk metal content (Thorngren et al. 2016). While atmospheric probes of large ( ≳ 2 . 4 R ⊕ ) Neptunes and sub-Neptunes revealed low atmospheric metallicities (Benneke et al. 2019a,b; Madhusudhan et al. 2023) in line with the expectations for a H 2 /He-dominated atmosphere, more diversity is expected in the direct vicinity of the radius valley, as H 2 /He atmospheres on smaller sub-Neptunes have greater susceptibility to atmospheric mass loss. We sought to explore this diversity of small  \nsub-Neptune atmospheres through a large JWST Cycle 2 program (GO 4098), by obtaining high-signal-to-noise ratio (S/N) transmission spectra of small, low-density planets with well-characterized masses and radii (Figure 1) across the 0.6 -5.5 µ m range with JWST's Single Object Slitless Spectroscopy (SOSS) mode of the Near Infrared Imager and Slitless Spectrograph (NIRISS) instrument (NIRISS/SOSS) and with the G395H grism of the Near Infrared Spectrograph (NIRSpec/G395H). Because of the small sizes and low equilibrium temperatures of the most promising candidates, transmission spectroscopy is the method of choice to characterize their atmospheric compositions.", '1.4. Introduction to the GJ 9827 system': 'GJ 9827 d is a 1.98 R ⊕ planet transiting a moderately active nearby (30 pc), bright (J magnitude of 7 . 984 ± 0 . 020) K7V dwarf (Prieto-Arranz et al. 2018). It is one of three planets in a 1:3:5 near-resonant system including two other closer-in super-Earths, GJ 9827 b (1.57 R ⊕ ) and GJ 9827 c (1.24 R ⊕ ), detected by the K2 survey (Rodriguez et al. 2018; Niraula et al. 2017). The coevolution of these three planets spanning the radius valley under the irradiation of the same host star makes GJ 9827 an ideal system to test our understanding of how atmosphere escape shapes planetary atmospheres and compositions.  \nRadial velocity observations with the Magellan II Planet Finder Spectrograph (Teske et al. 2018), the FIES, HARPS, and HARPS-N spectrographs (PrietoArranz et al. 2018; Rice et al. 2019), supported by later ESPRESSO observations (Passegger et al. 2024) revealed that while GJ 9827 b and c have rocky densities, GJ 9827 d has a lower density consistent with a volatile envelope. Indeed, even when the joint evolution of all three planets is taken into account, the absence of an atmosphere on GJ 9827 b and c does not preclude the presence of an atmosphere on GJ 9827 d (Owen & Campos Estrada 2020; Kosiarek et al. 2021; Bonomo et al. 2023).  \nThe mass and radius of GJ 9827d are consistent with compositions ranging from a fraction of a percent H 2 /He (Kosiarek et al. 2021; Passegger et al. 2024), or a mass fraction of HMMW volatiles of 5 to 30% with a mixedsupercritical envelope state for this T eq , A B =0 . 3 ≃ 620 K planet (Aguichine et al. 2021; Innes et al. 2023; Benneke et al. 2024). Atmospheric escape considerations, however, disfavor a hydrogen-dominated atmosphere, as the intense mass loss predicted in this case from both photoevaporation (Carleo et al. 2021), and corepowered mass loss (Gupta & Schlichting 2021), was faced with nondetections of Lyα , H α and He I using  \nhigh-resolution space-based and ground-based observations (Kasper et al. 2020; Carleo et al. 2021; Krishnamurthy et al. 2023).  \nWhile clouds often plague transmission spectra of subNeptune-size planets (Kreidberg et al. 2014; Wallack et al. 2024), the recently published HST /WFC3 transmission spectrum of GJ 9827 d exhibits an absorption feature at 1.4 µ m (Roy et al. 2023) which suggests that even if clouds partially mute the transmission spectrum, molecular features can still be detected. The analysis of the HST spectrum favored water over methane (which has a coinciding absorption band at 1.4 µ m, see Benneke et al. 2019b; Madhusudhan et al. 2023) as the source of the absorption feature, mainly from the shape of the wings of the molecular band. However, a broader transmission spectrum that would probe multiple rotational/vibrational bands of both molecules could disentangle their relative contributions to the observed absorption features in transmission (Benneke et al. 2024). Photospheric heterogeneities in the K7-type host GJ 9827 are not expected to create significant contamination via the transit light source effect (Rackham et al. 2018), lifting the ambiguity of the atmospheric or stellar origin of potential water absorption features (Roy et al. 2023).  \nHere we present the first JWST look at the atmosphere of GJ 9827 d using two NIRISS/SOSS transits. These observations are part of a large program that will additionally provide two transit observations with NIRSpec/G395H over the coming year. NIRISS/SOSS is sensitive to H 2 O and CH 4 absorption, with multiple absorption bands of each molecule across the 0.6 -2.8 µ m range. Furthermore, the resolution of NIRISS/SOSS opens up the possibility of constraining potential absorption from metastable He at 1083 nm (Fu et al. 2022).  \nOur study is outlined as follows. We present the NIRISS/SOSS observations and data reduction in Section 2, followed by the white- and spectroscopic lightcurve fitting procedure in Section 3. Our modeling of the atmospheric and stellar contributions to the spectrum is detailed in Section 4. We describe our atmospheric escape modeling in Section 5, and our planet structure modeling in Section 6. We present our results in Section 7, discuss them in Section 8, and summarize our conclusions in Section 9.', '2. OBSERVATIONS AND DATA REDUCTION': "In this study, we analyze two NIRISS/SOSS transit observations of GJ 9827 d using two independent data reduction pipelines. Throughout this study, we treat the supreme-SPOON reduction as our baseline and check that our results are not affected by the choice of reduc-  \nFigure 1. Mass-radius diagram showing the small planets in the sample of JWST Cycle 2 Program GO 4098 (triangle markers) and how they compare to the larger K218 b. We also display planets from the NASA Exoplanet Archive with 3 σ or better mass measurements, colored by their equilibrium temperatures for an Earth-like albedo of 0.3. The planets span both sides of the radius valley (gray). We highlight the compositional degeneracy in mass-radius space with constant composition lines calculated using smint for planets with H 2 /He envelopes or H 2 O in the supercritical+vapor state (mass-radius relationships from Lopez & Fortney (2014); Aguichine et al. (2021)). We display massradius curves for water-rich compositions of 10% and 50% water in supercritical+vapor state assuming an irradiation temperature of 670 K (representative of GJ 9827 d; blue solid lines), and isocomposition curves for 0.2% and 1%-by-mass solar-composition H 2 -dominated gas for a 1 Gyr planet receiving like GJ 9827 d approximately 35 times the Earth's insolation (navy dashed lines), both atop Earth-like cores. The positions of Venus and the Earth are indicated with black letters, and the light gray shaded area in the bottom right corner highlights planets smaller than expected from collisional stripping of planetary mantles (Marcus et al. 2010). With a mass of 3.02 M ⊕ and a radius of 1.98 R ⊕ , GJ 9827 d hugs the upper size end of the radius valley, in line with its high predicted hydrogen escape rates.  \n<!-- image -->  \ntion pipeline by comparing our transmission spectrum with those obtained from an independent NAMELESS reduction.", '2.1. Summary of the observations': "Two transits of GJ 9827d were observed with JWST's SOSS mode (Albert et al. 2023) of the NIRISS instrument (Doyon et al. 2023) as part of the Cycle 2 General Observer program GO 4098 (PI: Benneke; Co-PI: Evans-Soma). Each NIRISS/SOSS time series observation (TSO) lasted 3.4 hours. We chose the SUBSTRIP256 subarray which enables us to extract the stellar spectrum from both order 1 (0.8 to 2.85 µ m, average  \nresolving power R = 700) and order 2 (0.6 to 1.4 µ m, average resolving power R = 1400). We therefore probe the full 0.6 to 2.85 µ m range at medium resolution (Albert et al. 2023). Our target is bright in the J band, which makes NIRISS/SOSS the only JWST mode that can be leveraged to probe potential near-infrared water bands. We only observe two groups up the ramp to avoid saturation, with a total of 756 integrations over each TSO.  \nThe first transit (Visit 1) was observed on UTC November 10, 2023, and covered 0.4 hours of baseline before the transit, 1.3 hours in transit, and the following 1.7 hours out of transit. The second transit observation (Visit 2) occurred on UTC November 16, 2023, and captured 1.4 hours of pretransit baseline, the 1.3-hour transit, as well as 0.7 hours of posttransit baseline.", '2.2. supreme-SPOON reduction': "We reduce both SOSS visits using the supreme-SPOON pipeline 1 (Radica 2024; Feinstein et al. 2023; Radica et al. 2023), closely following the procedure in Radica et al. (2024) and Benneke et al. (2024). For each visit, we perform the standard supreme-SPOON Stage 1 and Stage 2 calibrations, including a group-level correction of column-correlated 1/ f noise, and a timedomain detection of cosmic rays specifically developed for ngroup=2 TSOs (Radica et al. 2024). For the background, we perform a piecewise correction (e.g., Lim et al. 2023; Fournier-Tondreau et al. 2023; Benneke et al. 2024), independently scaling the standard STScI SOSS background model on either side of the background 'step' to match the flux level of the observations. Moreover, due to the presence of a dispersed contaminant located above the third diffraction order of the target spectrum in both observations, we use pixels x ∈ [350 , 600] and y ∈ [230 , 250] to estimate the pre-step background level instead of the default region. We extract the 2D stellar spectra using a simple aperture extraction with a full width of 40 pixels around the trace since the order self-contamination is predicted to be negligible (Darveau-Bernier et al. 2022; Radica et al. 2022).", '2.3. NAMELESS reduction': 'We reduce the NIRISS/SOSS observations of GJ 9827 d using the NAMELESS pipeline (Feinstein et al. 2023; Coulombe et al. 2023) following the methods described in Coulombe et al. 2024 (submitted) and Benneke et al. (2024). Our reduction starts from the uncalibrated raw data and goes through the super-bias sub-  \ntraction, reference pixel correction, nonlinearity correction, ramp-fitting, and flat-fielding steps of the STScI jwst pipeline v1.12.5 (Bushouse et al. 2023). The jump detection step was automatically skipped by the pipeline as the integrations consist of only two groups. We then proceed with bad pixel correction where we flag and correct for pixels that show null, negative, or abnormally high counts at all integrations. The non-uniform background is subsequently subtracted by scaling independently the two regions of the STScI model background 2 which are separated by a sharp jump in flux around spectral pixel x ∼ 700, following the procedure presented in Lim et al. (2023). We also correct for cosmic rays by computing the running median of the flux in each pixel, using a window size of 10 integrations, and clipping any 4 σ outliers. Finally, we correct for the 1/f noise by scaling all columns of the first and second spectral orders independently and considering only pixels less than 30 pixels away from the center of the trace, as described in detail in Coulombe et al. 2024 (submitted). The light-curves are then extracted from the first and second spectral orders using a simple box aperture with a width of 36 pixels.', '3. EXOTEP LIGHT-CURVE FITTING': "We perform individual light-curve fits to each of the visits using the ExoTEP pipeline (Benneke et al. 2017, 2019a). We fit jointly the broadband light-curve fits of order 1 and order 2, to infer the system parameters, and then fix the planet's scaled semi-major axis a/R ⋆ and its impact parameter b to the best-fitting values from the joint white-light-curve fit when fitting wavelength-dependent light-curves. For NIRISS/SOSS, the supreme-SPOON extracted time-dependent spectra are used for the main analysis, and the resulting transmission spectrum is consistent with the one obtained with NAMELESS .", '3.1. Broadband Light-curve fitting': "For both TSOs, the NIRISS/SOSS white-light-curves for orders 1 and 2 exhibit the beating pattern that has been previously reported and is thought to correspond with the 2-4-minute cycling in the panel heating of the instrument electronics compartment (McElwain et al. 2023). These short-term variations are suppressed when detrending against the eigenvalues of the first two principal components extracted from the time series of detector images, which track changes in the trace morphology roughly corresponding to variations in the position (PC1) and width (PC2) of the trace as demonstrated in  \nthe eigenimages (see Coulombe et al. 2023). We found that even when detrending against the principal component analysis (PCA) components, residual correlated noise on longer timescales (10 to 30 min) is detected in the light-curves, which we tentatively attribute to stellar variability in the absence of corresponding variations in the recorded shifts of the spectrum or the trace width or position. We note, however, that since we are using a small number of groups up the ramp to accommodate the star's brightness, detector effects such as imperfect nonlinearity correction or outlier flagging that could affect time series observations of brighter targets might also have an impact on the extracted light-curves. We model this correlated noise source using the Gaussian Process (GP) kernel available in the celerite module that most closely reproduces a squared-exponential behavior (correlation between neighboring points over a characteristic timescale), the Matern 3/2 kernel.  \nOur final model combines the astrophysical (transit) component, as well as a systematics model. We use ExoTiC-LD (Grant & Wakeford 2022; Husser et al. 2013) to compute quadratic limb-darkening parameters, which we keep fixed (two parameters for order 1, two parameters for order 2) instead of fitting them, as the high impact parameter of GJ 9827 d poses a challenge to the joint fitting of the orbital and limb-darkening parameters. We use a custom SOSS order 2 throughput (O. Lim, priv. comm.) that covers all the extracted wavelengths, rather than the default ExoTiC-LD instrument throughput. Our systematics model has six free parameters: a linear trend with time, which is frequently observed in other JWST datasets (e.g., Coulombe et al. 2023; Feinstein et al. 2023), a linear function of the first two PCA eigenvalues approximately tracing vertical shifts and variations of the width of the trace, and residual correlated noise captured by the Matern 3/2 GP term, where we fit the correlation strength σ GP and timescale τ GP . We fit the white-light-curves of orders 1 and 2 together in order to extract the maximum information content on the orbital parameters of the planet that are fitted: the impact parameter b , scaled semimajor axis a/R ⋆ and the transit time T c .  \nBeyond the previously mentioned systematics, visual inspection of the first NIRISS/SOSS visit reveals a jump in the light-curve around BJD TDB = 2460258 . 91 (in the transit) with a sudden increase followed by a decrease in flux, and a sudden decrease in flux close to BJD TDB = 2460258 . 955 (in the post-transit baseline). Besides, the light-curve of the second NIRISS/SOSS visit displays a short in-transit bump, which we investigate in more detail and attribute to a stellar flare (see Section C). To model these visit-specific effects independently, we  \nTable 1. Stellar and planetary parameters adopted or refined in this work. The ephemeris is obtained by fitting a constant period and transit time to the observed Kepler , HST , Spitzer , and NIRISS/SOSS transit times together. The reference transit time is anchored to the epoch of the last JWST NIRISS/SOSS transit. The stellar luminosity is calculated from the literature values of the effective temperature and radius of GJ 9827. The planetary orbital separation, insolation, and equilibrium temperature values are derived from the scaled semi-major axis constraints obtained with NIRISS/SOSS, and using the literature values of the stellar radius and effective temperature. The planetary radius is obtained from the planet-to-star radius ratio fitted on the NIRISS/SOSS order 1 white light curves combined with the stellar radius constraint.  \nperform separate white-light-curve fits for the first and second visits.  \nWhile the small in-transit bump observed in visit 2 does not bias the determination of the system parameters from the light-curve as the egress remains unaffected by the flare, the feature in the transit for visit 1 is much broader and hinders independent determination of b and a/R ⋆ , as it cannot be efficiently masked while obtaining sufficient phase coverage of the transit. Therefore, we first fit the cleaner light-curve of the second visit and extract from their posterior distributions the covariance matrix for b and a/R ⋆ , which we use as a correlated Gaussian prior on these two parameters when fitting the light-curve of the first visit (Figure 2).  \nFor the second visit, we model the in-transit bump by adding a narrow Gaussian bump to the systematics model, which has three more parameters: the center of the Gaussian t sc [BJD TDB ], its width σ sc and its amplitude A sc (see Roy et al. 2024, in prep).  \nBefore performing the final fit to the white light curves of visit 1, we obtained a Gaussian prior on its transit time T c , 1 from a simple fit to the portions of the light-curve in both orders that appeared least affected by stellar or instrumental nuisance signals (excluding BJD ∈ [2460258 . 85 , 2460258 . 87] and BJD ≥ 2460258 . 955). For this simple fit, the systematics model did not include a GP and only accounted for a linear detrending with time and with the eigenvalues from the first two principal components. Finally, we ran a full fit to the  \nFigure 2. Top panels: Results from the white light curve fit of the supreme-SPOON reduction for Visits 1 (left) and 2 (right), and fitted systematics models. The order 1 systematics-corrected light-curves are shown (teal points), binned for visual purposes, along with the best-fit transit model (black) and the systematics model including the GP component (red, with the GP scaled for visualization purposes). The vertical blue shaded regions indicate the three candidate flares with the highest signal-to-noise spectra and following typical flare evolution profiles (see Section C). Bottom panels: Spectroscopic light-curve fits to the two NIRISS/SOSS transits of GJ 9827 d from the supreme-SPOON reduction. The left panel corresponds to the first visit and the right panel to the second visit. We show the detrended broadband light-curve for order 2, and seven spectroscopic light-curves from order 1, from top to bottom. The best-fit astrophysical model is shown for each light-curve (black) and light-curves are offset relative to each other for clarity.  \n<!-- image -->  \nentire light-curve, introducing again the Matern 3/2 GP component, and using a covariance matrix to encode our prior on the values of T c , 1 , b , and a/R ⋆ (Figure 2).  \nWe find consistent planetary parameters across both visits (Table A1) and in agreement with literature values (Table 1; Section A). We also obtain an updated ephemeris for GJ 9827 d (see Table 1; Section A), and do not detect significant transit-timing variations (TTVs) using the two additional transit time measurements (see Figure A1; Section A).", '3.2. Spectroscopic light-curve fitting': 'We obtain the transmission spectrum of GJ 9827 d for each visit by fitting the spectroscopic light-curves using ExoTEP (Benneke et al. 2017), and fixing the orbital parameters ( b , a/R ⋆ ), and the mid-transit time T c to the best-fit parameters from the white light curve. We bin the extracted 2D spectra to obtain three versions of the transmission spectrum. The first version of the spectrum, at a resolving power of 100, is used in the retrieval analyses as it provides a good compro-  \nse between computationally expensive retrievals on a higher-resolution spectrum and loss of spectral information for heavier wavelength binning (Welbanks et al., in prep). We also produce two higher-resolution spectra where we only extracted the order 1 spectrum around the metastable He I triplet to search for any He absorption feature. For these He search analyses, we use one version of the spectrum obtained by fitting the transit depth in every detector column and one where we bin together two neighboring columns. We verify that all three versions of the spectrum agree in the spectral region around the metastable He triplet where the highresolution spectra were obtained (Figure A2).  \nSimilarly to the white light curve analysis, we fit jointly the astrophysical model described by the planetto-star radius ratio R p /R ⋆ , and a systematics model including parameters that capture instrumental and astrophysical noise, to each spectroscopic light-curve. Because of the high impact parameter of the planet (Table 1), we keep the limb-darkening parameters fixed to the 1D quadratic parameters computed in each bin by  \nExoTIC-LD using the stellar parameters of GJ 9827. We fit to each light-curve a scatter term s , the two parameters of a slope with time c + vt , and slopes with the eigenvalues of the first two principal components, v PC1 and v PC2 . For visit 2, we additionally fit the amplitude and width of a Gaussian describing the bump, with the feature center time fixed to the best-fit value from the white light curve fit.  \nWe try fitting GP parameters to each spectroscopic light-curve but find that they are unconstrained in most spectroscopic bins and lead to overestimated transit depth uncertainties. Rather, we elect to use the recorded GPmodel hyperparameters and the corresponding mean GP model of the residuals from the best fit to the white light curve. For the first visit, we fix the Matern 3/2 GP lengthscale to the best-fit white light curve values and model the residuals in each spectroscopic bin by only fitting for the correlation strength parameter σ over a narrow range of values (factor of 0 to 10) around the best-fitting value from the white light curve fit. For the second visit, where the structure of the residuals in each spectroscopic bin match more closely those of the white light curve, we use the mean GP model computed from the GP parameters of the white light curve best-fit systematics model and scale it by fitting a multiplicative factor in each bin. We check that the residuals from the fitted spectroscopic light-curves (Figure 2) follow the distribution expected from white noise using Allan plots.  \nThe visit 2 spectrum is robust to the choice of GP model for the spectroscopic light-curve fit as the scaled white light curve model or the scaled best-fit GP model parameters to the white light curve, and we report the spectrum obtained using the less-flexible approach (e.g., Radica et al. 2024). For both visits, we confirm that the spectrum we obtain with this approach is consistent within less than 1 σ with the spectrum obtained from the more conservative approach of excluding the regions of the light-curve that are directly within the in-transit bump features.  \nWe repeat the same white- and spectroscopic lightcurve fitting procedure for the supreme-SPOON reduction of both visits, which produces a spectrum consistent with NAMELESS at better than 1 σ (Figure A3). We note, however, that there are slight differences between the two reductions in how well the extracted principal components correlate to the white light curve residuals. The correlation is stronger for the NAMELESS reductions, which leads to a lower residual scatter. The overall slope in the white light curve baseline is also different between the two reductions, but it does not affect the resulting spectrum as the small differences (with a mean differ-  \nnce of 11 ppm overall) are absorbed by the slope and GP parameters.  \nOur spectrum is also in excellent agreement with the spectrum obtained with HST (Figure A4) from combining 10 transits of GJ 9827 d (Roy et al. 2023), with a small vertical shift which we attribute to the fact that the orbital parameters b and a/R ⋆ were fixed to literature values (Niraula et al. 2017) due to partial ingress/egress coverage when fitting the HST data.  \nWe calculate the final spectrum by averaging the spectra from the two visits and use the spectrum from the supreme-SPOON reduction for the atmospheric analysis (Table A2).', '3.3. Interpretation of the observed stellar variability features': "We study the wavelength dependence and time evolution of the residual features in the light-curves that are captured by the GP in our fit (Section C). Although the time resolution offered by NIRISS/SOSS does not allow us to probe the detailed time evolution of each feature (Figure A5), we prefer flares to stellar spot crossings because of the 'fast-rise, exponential decay' (FRED) evolution observed in the time series of flare effective temperatures that we derive from the observations for three of the candidate flare events, labeled F1, F3, and F4 (Figures 2 and A5). We do not detect H α variability, which was found to be a telltale flare diagnostic for the much cooler star TRAPPIST-1 (Lim et al. 2023; Howard et al. 2023), and attribute this to the difference in brightness and spectral energy distribution between the M7.5V star TRAPPIST-1 and the K7V star GJ 9827, where flares may be characterized by continuum enhancements rather than H α variations (see Figure A6; discussion in Section C). The short-wavelength coverage of NIRISS/SOSS enables us to extract, for each of the candidate flare events, the flare energy in the TESS bandpass, and to leverage relations between the TESS flare energy and the near-UV and far-UV (NUV, FUV) energies to estimate the NUV and FUV outputs corresponding to each candidate flare (Table A3).", '4. ATMOSPHERE AND STELLAR CONTAMINATION MODELING': 'We use SCARLET and POSEIDON to perform 1D retrievals of the atmospheric makeup of GJ 9827 d from the NIRISS/SOSS and WFC3 transmission spectra. We also consider the potential impact of unocculted stellar heterogeneities on the inferred atmospheric properties by fitting the spectrum assuming that it can be entirely explained by stellar contamination, or accounting for the contributions of both the planetary atmosphere and any  \npotential stellar surface heterogeneities on the observed features.', '4.1. Atmosphere modeling with SCARLET': 'We perform both free chemistry and chemically consistent retrievals over the spectrum of GJ 9827 d using the SCARLET atmospheric retrieval framework (Benneke & Seager 2012, 2013; Benneke 2015; Benneke et al. 2019a,b; Pelletier et al. 2021; Piaulet et al. 2023).  \nThe version of SCARLET used in this study leverages several upgrades over previously published work, including the implementation of the joint fitting of a stellar contamination contribution to the spectrum using the MSG module (Townsend & Lopez 2023) alongside the atmosphere parameters, and the introduction of partial aerosol coverage using a cloud fraction (Welbanks & Madhusudhan 2019). In chemically consistent retrievals, we also run cases where we allow the CH 4 and/or NH 3 abundance to be fitted freely independent of the chemical equilibrium prediction, to test the robustness of our results and potential disequilibrium depletions and avoid unphysically low C/O ratios if e.g. methane is not detected.  \nWe explored a range of assumptions to explain the spectrum of GJ 9827 d, from free retrievals with wellmixed abundances of the fitted species to chemical equilibrium retrievals that follow the predicted molecular abundances in thermochemical equilibrium for the fitted isothermal temperature, atmospheric metallicity, and C/O ratio, with or without a freely fitted CH 4 or NH 3 abundance to capture potential disequilibrium depletion.  \nFor free retrievals, the forward model of SCARLET performs the radiative transfer and iterative hydrostatic equilibrium calculation for the atmospheric composition and temperature structure prescribed by the sampled parameters, while for chemically consistent retrievals the molecular abundances are set not only by the sampled metallicity and C/O ratio, but also by the local pressure and temperature conditions in each layer of the atmosphere. Our free retrieval includes H 2 , He, N 2 , HCN, H 2 O, CO, CO 2 , CH 4 , NH 3 , H 2 S, and for absorbing species, we use the HELIOS-K computed cross sections of H 2 O (Polyansky et al. 2018), CO (Hargreaves et al. 2019), CO 2 (Yurchenko et al. 2020), CH 4 (Hargreaves et al. 2020), HCN (Harris et al. 2006), NH 3 (Coles et al. 2019), H 2 S (Azzam et al. 2016).  \nFor each combination of model parameters, the radius at the reference pressure of 10 mbar is fitted to find the value that provides the best match between the model spectrum and the observations, with the hydro-  \nstatic equilibrium and radiative transfer steps iterated over until the optimal radius value is found.  \nThe chemical composition of the planet is parameterized either by the fitted abundances (free retrieval) or by the log atmospheric metallicity, C/O ratio, and CH 4 abundances (chemically consistent retrievals). We fit for the temperature T p of the photosphere which is sampled by our observations. For the aerosols, we use either a gray cloud and hazes or a more complex Mie cloud parameterization (Benneke et al. 2019a). In the first scenario, clouds and hazes are parameterized respectively by the log of the cloud-top pressure (approximately τ = 1 level), and c haze , the slope enhancement parameter relative to Rayleigh scattering, caused by small particles. In the Mie cloud scenario, we fit for the scale height of the cloud relative to the atmosphere, H rel , for the density of cloud particles n and for the mean particle size µ . We allow for patchy clouds with the f cloud parameters, where clouds are patchy if f cloud < 1. The parameter space is sampled using multi-ellipsoid nested sampling as implemented in the nestle module (Skilling 2004, 2006). Models are computed at a resolving power of 15,625 and then convolved to the resolution of the observed spectrum for the likelihood calculation. We additionally fit for a potential offset between the Order 1 and Order 2 transmission spectrum and between the Order 1 NIRISS/SOSS and HST /WFC3 spectra.  \nFor the temperature profile, we fit the temperature of the photosphere of the atmosphere at the terminator, which is representative of the region probed by our transmission observations. In free retrievals, we fit for the abundances of the molecular species that have spectral features within the NIRISS/SOSS bandpass, with wide log-uniform priors on the volume mixing ratios (lower bound at 10 -10 ): H 2 O, CH 4 , and NH 3 , as well as CO, CO 2 , H 2 S, and HCN which have much weaker expected opacity over the NIRISS/SOSS wavelength range (few absorption bands and/or low expected volume mixing ratios). We ran retrievals where the background, or filler gas is H 2 /He (enforcing a Jupiter-like He/H 2 = 0.157), H 2 O, or N 2 .', '4.2. Atmosphere modeling with POSEIDON': 'We use the combined transmission spectrum of GJ 9827 d using both JWST/NIRISS SOSS data and HST/WFC3 data from Roy et al. (2023) to infer properties of the atmosphere. We do this by also conducting atmospheric retrievals with the open-source retrieval code POSEIDON (MacDonald & Madhusudhan 2017; MacDonald 2023), which utilizes the nested sampling algorithm PyMultiNest (Feroz et al. 2009; Buchner et al. 2014) to investigate the potential atmospheric  \nmodels, from which we use 2000 live points for each of our models. The descriptions of the radiative transfer technique, forward atmospheric modeling, and the opacity database that POSEIDON employs can be found in MacDonald & Lewis (2022). From our data, we produce spectra from 0.6 to 2.9 µ m at R = 20000, which is convolved with a Gaussian kernel to the native resolution of NIRISS/SOSS (R ≈ 650) and HST/WFC3 (R ≈ 150), respectively, and then multiplied by the instrument sensitivity function and subsequently binned down to the wavelength spacing of our data.  \nWe consider a wide range of atmospheric models, including a flat line model, pure H 2 O atmospheres, multigas atmospheres, and with or without stellar contamination (see Section 4.5). Our multigas models include the following gases: H 2 , He, N 2 , HCN, H 2 O, CO, CO 2 , CH 4 , NH 3 , and H 2 S. We note that our POSEIDON multigas models also use a permutation-invariant centered-logratio (CLR) prior (Benneke & Seager 2012) for the volume mixing ratio, meaning any molecule is theoretically equally likely to be the dominant gas. We fit for a reference radius, R p , ref , at a reference pressure of 10 bars, assuming an isothermal atmospheric temperature. All of our non-flat models include both clouds and hazes, where we assume an optically-thick gray opacity, and all layers deeper than P surf are set to infinite opacity. We also use a two-parameter power-law prescription for hazes (MacDonald & Madhusudhan 2017) in all of our multigas atmospheric models. The priors we use in our retrievals are as follows: R p , ref (R ⊕ ) = U (0.60, 1.15), T (K) = U (100, 1000), log a haze = U (-4, 8), γ haze = U (-20, 2), log P cloud = U (-7, 2), log X = CLR (-12, 0).', '4.3. Stellar-contamination-only retrievals with stctm': "We use the open-source stctm module (PiauletGhorayeb 2024) previously applied to TRAPPIST-1 b (Lim et al. 2023), TRAPPIST-1 g (Benneke et al., submitted), and GJ 9827 d (Roy et al. 2023) to model the case where the planet's atmospheric signature would have no wavelength dependence. This corresponds to a scenario where the planetary atmosphere contribution to GJ 9827 d's transmission spectrum would be featureless, e.g., from high-altitude opaque uniform clouds, or in a bare rock scenario which is unlikely given the planet's low bulk density.  \nIn order to check whether stellar surface heterogeneities can completely explain the water features in GJ 9827 d's spectrum, we use stctm to perform stellar contamination-only retrievals of the transmission spectrum assuming that (1) only spots contribute, and (2) that both spots and faculae play a role (parameterization described in Appendix Section D).  \nIn the spots-only case, we fit for ∆ T spot < 0 (see Appendix Section D), for the spot covering fraction f spot , and the planet's wavelength-independent transit depth D . In the spots+faculae case, we fit for two additional parameters: ∆ T fac > 0 and f fac which describe the temperature and covering fraction of faculae. To obtain the specific intensity I λ of the spots, faculae, and photosphere, we use the MSG module (Townsend & Lopez 2023) which interpolates across the PHOENIX grid of stellar atmosphere models (Husser et al. 2013). To speed up each likelihood evaluation, we pre-compute a grid of stellar models at a resolving power of 10,000 with MSG , which match the star's surface gravity log g , and span the effective temperatures T eff in small intervals of 20 K. In the retrieval, we use the stellar model within the grid corresponding to the closest T eff to the sampled value of T spot , T fac or T phot .  \nThe retrievals leverage the emcee implementation of Markov-Chain Monte Carlo sampling to explore the parameter space (Foreman-Mackey et al. 2013). We use 20 times as many walkers as there are fitted parameters, and run the chains for 5,000 steps, discarding 60% as burn-in prior to obtaining the final posterior distributions on the heterogeneity parameters. We use uniform priors on the heterogeneity fractions and ∆ T het , with | ∆ T het | > 100 as otherwise the heterogeneity spectra are so similar to the photosphere spectrum that any heterogeneity fractions can be invoked. We jointly fit for the photosphere temperature T phot to propagate its uncertainty, with a Gaussian prior informed by the literature constraint (Table 1, Passegger et al. 2024). The faculae fraction and temperature being largely unconstrained in the second retrieval configuration, without improving the quality of the fit to the data, we only present results for the spot-only fit.", '4.4. Joint stellar contamination and planetary atmosphere retrievals with SCARLET': 'We perform atmospheric retrievals with SCARLET using a conservative treatment of the potential transit light source (TLS) effect on the spectrum where we marginalize over any contributions consistent with the observations. We use the SCARLET implementation of joint atmosphere and stellar contamination retrievals that was developed for the study of HAT-P-18 b (Fournier-Tondreau et al. 2023), and update the stellar modeling query component to use the MSG module instead of the standard pysynphot implementation.  \nFor these joint retrievals, we therefore simultaneously sample the parameters of the planetary atmosphere and of the stellar surface, as parameterized in the stctm implementation. The stellar contamination factor ϵ λ, het is  \nmultiplied by the predicted transmission spectrum given the sampled atmosphere parameters prior to the likelihood evaluation. We test cases where we consider no stellar contamination, spots only (with a single population having a shared T spot ), and both spots and faculae with an associated temperature for each heterogeneity population.', '4.5. Joint stellar contamination and planetary atmosphere retrievals with POSEIDON': "We further use POSEIDON to determine the degree to which the TLS effect may be able to explain GJ 9827 d's transmission spectrum in the absence of an atmosphere. POSEIDON produces model spectra by multiplying a bare-rock transmission term, (R p /R ∗ ) 2 , by the wavelength-dependent stellar contamination factor from two distinct stellar heterogeneities (see Equation D1). Using the PyMSG package, POSEIDON calculates the spectra of the active regions by interpolating across the PHOENIX grid of stellar atmosphere models (Husser et al. 2013). We add five additional free parameters to each of our TLS models, the priors of which are as follows: f fac = U (0.0, 0.5), f spot = U (0.0, 0.5), T fac (K) = U (T ∗ -36, 1.2T ∗ ), T spot (K) = U (2300, T ∗ +36), T phot (K) = U (T ∗ , 12), where T ∗ = 4236 ± 12. The surface gravity of the stellar regions is fixed to log g = 4.719 ± 0.02 (Passegger et al. 2024).", '5. ATMOSPHERIC ESCAPE MODELING': 'We examine hydrodynamic escape in more detail, using the CHAIN (Cloudy e Hydro Ancora INsieme) 1D hydrodynamic model introduced in Kubyshkina et al. (2024) (see details of the modeling setup in Appendix Section E).  \nWe compute mass loss rates, as well as temperature, pressure, density, and ionization profiles with height for several composition cases. We take as the reference case a H 2 /He-dominated atmosphere with a metallicity matching that of the star and also run models for water-enriched H 2 /He atmospheres (with varying mass fractions of water mixed with the hydrogen and helium, see Egger et al. 2024 for more detail), or H 2 /He atmospheres with increased metallicity relative to solar (following the same methodology as Linssen et al. 2024). The water- and metal-enriched simulations allow us to account for the effect of enrichment in heavier species on the stability of the atmosphere against escape.', '6. PLANETARY STRUCTURE MODELING': "We use smint (Piaulet et al. 2021, 2023), an opensource Bayesian estimator of planetary compositions, to constrain the properties of GJ 9827 d in the scenario  \nwhere water is the only volatile present on top of a rocky core with an unconstrained ratio of silicates to iron. The model used accounts for the various phase transitions of water depending on the temperature and pressure conditions (Aguichine et al. 2021) which, for GJ 9827 d, corresponds to water in the gas, then in the supercritical phase from low to high pressures due to its irradiation level.  \nWe use a Gaussian prior on the irradiation temperature of GJ 9827 d ( T irr = 674 ± 29 K), calculated from Monte Carlo sampling of the stellar effective temperature, radius, and the planetary orbital distance (Table 1) and on the planet mass (Passegger et al. 2024). We uniformly sample the fraction of iron in the iron core + silicate mantle interior of the planet, f ' core ( f ' core = 0 in the absence of an iron core, f ' core = 325 for an Earthlike interior composition) and use a uniform prior on the water mass fraction f H 2 O . The parameter space is sampled using emcee (Foreman-Mackey et al. 2013). We run the chains for 2,000 steps and discard the first 60% as burn-in to obtain the final posterior distributions.", '7. RESULTS': "We perform retrievals to jointly interpret the HST/WFC3 (Roy et al. 2023) and our NIRISS/SOSS transmission spectrum of GJ 9827 d (Figure 3), which cannot be explained by stellar contamination alone (Figure A7). Atmospheric retrievals favor a high mean molecular weight, water-rich atmosphere (Figures A8,4). We also search for a signature of escaping metastable helium by leveraging the high spectral resolution in NIRISS/SOSS order 1 but do not detect any significant signal. Our hydrodynamic simulations predict atmosphere mass loss rates much higher than our measured upper limit. Overall, the data are best explained by a 'steam world' atmosphere with a high water abundance.", '7.1. Constraints on the atmospheric composition': "We fit the JWST/NIRISS SOSS + HST/WFC3 transmission spectra of GJ 9827d jointly, as the large number of transits observed with HST (10) yielded small enough error bars to be comparable in precision with the NIRISS/SOSS spectrum (Figure A4). We find that the absorption features in the transmission spectrum of GJ 9827 d can be robustly attributed to water vapor in its atmosphere, as stellar spots alone cannot reproduce the amplitude of the observed features (Figure A7), in agreement with the HST results (Roy et al. 2023). We do not detect any other molecular species. Most of the posterior probability mass of free retrievals regardless of the assumed background gas (Figure A8) and of  \nFigure 3. Results from the free retrieval performed on the JWST/NIRISS SOSS + HST/WFC3 spectrum of GJ 9827 d. Left panels: The semi-transparent lines are sample spectra from the posterior distributions of the SCARLET free retrievals when only the atmosphere (top panel) or both the atmosphere and stellar contamination component (bottom two panels) are accounted for. The black line is the best-fit model, and the black (gray) points are the measured NIRISS/SOSS (HST/WFC3) transit depths. The retrievals were performed using the R=100 version of the JWST/NIRISS SOSS spectrum, and a binned version is shown for visual comparison purposes. Even when stellar contamination is accounted for, the planetary atmosphere component of the model captures the H 2 O absorption features, and the constraints on the H 2 O abundance are not affected. The model including the contribution from spots only introduces an additional slope and does not affect the retrieval results. Right panels: Posterior distributions obtained on H 2 O in both retrieval setups, and on the spot covering fraction and temperature in the joint retrieval. For the POSEIDON retrievals (CLR prior) we only show in this figure the samples where the spot is at most 800 K cooler than the photosphere, for clarity (see full posteriors in Figure A8). Our inference of a water-rich, steam atmosphere is robust to the marginalization over stellar contamination contributions, and to the choice of priors (H 2 /He as the background gas in blue, CLR prior in pink).  \n<!-- image -->  \nthe chemically-consistent retrievals (Figure 4) lies in a metal-enriched regime with dozens of percent of the atmosphere mass contained in metals (Figure 4).  \nWe find that a high mean-molecular weight waterrich atmosphere can reproduce the spectrum even without any cloud opacity (Figure 4). Free retrievals performed with SCARLET and POSEIDON favor high water abundances, with tentative upper limits on the CH 4 and NH 3 abundances driven by the nondetection of the 2.3 µ m CH 4 features (other features overlap with water bands), and weak constraints on the 1.5 and 2.0 µ m ammonia features (Figure 3). CO and CO 2 abundances are not significantly constrained by our observations (Fig-  \nes A9, A10). Retrievals with CLR priors specifically favor H 2 O rather than H 2 as the dominant species in the atmosphere.  \nWith H 2 O as the only O-bearing molecular species with detectable bands in the NIRISS/SOSS bandpass and no detection of any C-bearing species, our constraints, especially on the C/H ratio, remain wide (Table 2), while the measured O/H ratio is precise even when marginalizing over undetected CO or CO 2 . Potential SO 2 in the atmosphere (with no significant absorption in the NIRISS/SOSS bandpass) may alter our O/H estimate. The most agnostic retrieval setup run (CLR free retrieval) measures a mean atmospheric molecular  \n<!-- image -->  \nFigure 4. Constraints on the atmospheric composition from the POSEIDON retrieval with the agnostic (CLR) prior on the background gas (pink), the SCARLET retrieval with H 2 /He assumed to be the background gas (blue), and SCARLET chemically-consistent retrieval (orange). Left panel: Contours representing the 0.5, 1, 1.5, and 2 σ levels in the joint posterior probability distribution on the atmospheric metallicity and the gray cloud top pressure, for all three retrievals. Right panel: Kernel density estimations of the atmospheric metal mass fraction Z atm obtained from all posterior samples, for each retrieval setup.  \n<!-- image -->  \nweight that closely matches that of water (18 . 01 +0 . 20 -0 . 61 ) and estimates that about four times as much mass is contained in O as in H.  \nOur constraints on the atmospheric composition are robust to the choice of filler gas in the SCARLET retrievals: H 2 /He (top panels in Figure A8), H 2 O, N 2 (not shown), or an equal prior on any molecule being the filler species (centered-log-ratios retrievals performed with POSEIDON, see bottom panels in Figure A8). The interpretation is consistent across the assumptions made, from chemical equilibrium to free retrievals, with small variations related to the prior on how molecular abundances vary with height (Figures A8 and 4). We note that even when using log uniform priors for the volume mixing ratio (H 2 /He assumed as the background 'filler' gas), we find similar strong evidence of a H 2 O-dominated atmosphere, in full agreement with the SCARLET results.  \nThe Mie cloud parameterization of particle size distributions (Benneke et al. 2019a) was found not to be warranted by the quality of our data, from a model comparison standpoint, and all the retrievals presented model aerosols as a gray cloud and a potential haze slope in all the results discussed. Our data only marginally ( ∼ 1 . 9 σ ) favor CH 4 depletion compared to chemical equilibrium (Table A4).", '7.2. Impact of stellar contamination and instrument offsets': "We investigate the impact of accounting for potential stellar surface heterogeneities on the results from our  \nretrievals. We show that for both free and chemically consistent retrievals, the addition of the contamination signal from spots (cooler regions on the stellar surface) and faculae (hotter plages) is not significantly favored from a Bayesian model comparison standpoint (Table A4). Further, we demonstrate that the addition of spots, or both spots and faculae, does not alter our inference of atmospheric properties (Figure A8). This finding is in line with previous HST /WFC3 results (Roy et al. 2023) which rule out stellar contamination as the source for the 1.4 µ m spectral band that was detected with HST .  \nIn particular, contrary to M dwarfs where spots can often easily reproduce planetary atmosphere-like water bands, the K7 dwarf GJ 9827 would require unrealistically large covering fractions of spots over ∼ 1000 K cooler than the quiet photosphere to produce the observed 1.4 µ m water absorption features (see modeling performed in Roy et al. 2023, and Figure A7). Modeling performed with POSEIDON for TLS-only models agrees with the stctm models and cannot reproduce the observed features.  \nWhen stellar contamination is considered in addition to a planetary atmosphere, testing for spots, spots+faculae, and without any stellar contamination for Bayesian comparison, we only find weak Bayesian evidence for the spot-only model, while we find no evidence of stellar faculae (Table A4). We try independently fitting each NIRISS/SOSS visit with or without spots or faculae to be robust to the presence of different stellar contamination contributions in each visit's spectrum, but the addition of stellar heterogeneity is  \nTable 2. Atmospheric retrieval results for GJ 9827 d. We quote results from the free retrievals including one population of stellar heterogeneity (spots) with POSEIDON (CLR prior case, for the samples where T spot is no more than 800 K cooler than T phot , star ) and SCARLET (H 2 /He as filler gas case). The uncertainties represent the 68% confidence intervals for all two-sided constraints and 95% (2 σ ) upper or lower limits otherwise.  \ndisfavored from a model comparison standpoint in each case.  \nWe also fit for a vertical offset between the NIRISS/SOSS order 1 and order 2 transmission spectra, and we find that the retrieved offset is very small compared to the error bars on the individual transit depths. For example, the H 2 /He background free retrieval with stellar spots finds an offset of -4 +21 -20 ppm. However, the addition of an offset (retrieved to be 15 . 7 +5 . 0 -4 . 6 for the atmosphere + spots case) is essential to jointly fit the JWST and HST data due to the difference between the assumed (HST) and fitted (JWST) orbital parameters of the planet (see Table 2).", '7.3. Search for escaping metastable helium': "The first order of NIRISS/SOSS covers the metastable He I triplet, which could show up in the transmission spectrum as a result of an extended escaping H/He atmosphere on GJ 9827 d. We sought to leverage the resolving power of NIRISS/SOSS to search for a He I signature. We produced two new spectral extractions, one at full resolution, and one at 2-pixel resolution, in the vicinity of the He I triplet (Figure 5). We fit the full resolution spectrum a Gaussian with a fixed width (0.75 ˚ A) convolved at the resolution of NIRISS/SOSS (700), and at the expected location of the He I triplet absorption (10833.33 ˚ A), to constrain the sensitivity of our observations to a potential signal.  \nFigure 5. Search for metastable He absorption around 10,833 ˚ A in both NIRISS/SOSS transits of GJ 9827 d. The top and bottom panels respectively show the spectrum fitted after extraction at pixel resolution, and 2-pixel resolution, for each of the two visits ( supreme-SPOON reduction). We do not detect significant absorption from the He I metastable triplet, in line with more sensitive ground-based high-resolution observations.  \n<!-- image -->  \nWe do not report any detection of the helium triplet in the NIRISS/SOSS full-resolution spectrum of GJ 9827 d. We obtain a constraint on the amplitude of the helium absorption signal in our data which corresponds to an upper limit on the amplitude of a helium absorption signal of 0 . 56 ± 0 . 23 % at high spectral resolution (80,000). The lack of metastable helium absorption on GJ 9827 d is surprising, as planets in similar irradiation conditions to GJ 9827 d around K-type stars are more likely to have helium particles in their metastable state leading to the detection of extended thermosphere and escaping material (Oklopˇci'c 2019). Still, our upper limit is consistent with more constraining nondetections of escaping metastable He and neutral H obtained from the ground with Keck/NIRPSEC, CARMENES or the InfraRed Doppler (IRD) instrument on the Subaru telescope (Kasper et al. 2020; Carleo et al. 2021; Krishnamurthy et al. 2023).", '7.4. Hydrodynamic escape predictions': 'For the hydrogen-helium dominated atmosphere simulation with a stellar metallicity, our model predicts an atmospheric mass loss rate of ∼ 2 × 10 10 g s -1 , and the escape is dominated by photoionization heating from atomic hydrogen. The ion fraction in the H 2 /He dominated atmosphere remains relatively low ( f ion is only ∼ 0 . 2 at the Roche lobe), suggesting that hydrogen escape should be detectable in Ly α line (see e.g., Kubyshkina et al. 2022).  \nFor water-enriched atmosphere cases compatible with the range of atmospheric water enrichments indicated by  \nthe transmission spectrum of GJ 9827 d (Figure 4), we obtain a reduction in the atmospheric mass loss rates by a factor of ∼ 2 (for a water mass fraction of 70%) and up to an order of magnitude (factor of ∼ 8 for a water mass fraction of 95%). The reduction is similar in the highmetallicity atmosphere cases, with slightly lower mass loss rates compared to water-rich atmospheres for the same atmospheric mean molecular weight, due to heating from metal ions. While the increase in mean molecular weight is the main driver for the decrease in mass loss rate, oxygen line cooling (for the water-rich cases) and cooling via the metal lines (for the high-metallicity cases) can contribute up to 90% of the total radiative cooling in metal-enriched atmospheres.  \nWe also find that even though a large fraction of the escaping atmospheric material still consists of atomic hydrogen in water-rich atmospheres, the increasing mean molecular weight leads to more compact atmospheres with narrower ionization fronts (the ion fraction increases faster with altitude). Therefore, neutral hydrogen particles are limited to low altitudes, and the total ion fraction in the outflow increases by up to a factor of a few hundred at the sonic point. We note, however, that our models that use the Cloudy line lists are not conclusive as to the impact of increased metallicity on the detectability of the metastable He feature beyond the fact that helium itself would be less abundant in atmospheres composed largely of water, or of metals more generally.  \nFinally, our models predict that for water-rich atmospheres (but much less so for high-metallicity atmospheres), a large fraction of the escaping material is represented by atomic oxygen supplied by the photodissociation of water molecules above µ bar pressure level. At high altitude, most of the oxygen is expected to be in ion form, with O + and O 2+ making up as much as 70% of the outflow by mass at the sonic point (up to ∼ 3 orders of magnitude more than in H 2 /He dominated atmospheres; see Egger et al. 2024). It remains unclear, however, whether these oxygen species would escape along with the hydrogen (see Section E).', '7.5. Structure modeling results': "If the entire volatile content of GJ 9827 d was contained within a pure water vapor/supercritical envelope (Figure A11), we find that its water mass fraction would be 32 ± 10% when marginalizing over the full range of plausible iron fractions in the planetary interior. This constrains to approximately 60-80% the total mass of the planet locked in its rock/iron core, assuming the mantle has completely solidified.  \nOur water mass fraction is consistent with the 5-30% reported by Aguichine et al. (2021) for GJ 9827 d. The differences between the two studies can be traced back to the fact that Aguichine et al. (2021) used a previous determination of the planetary parameters with a higher mass for GJ 9827 d (Rice et al. 2019), and did not marginalize over the uncertainty on the planet's irradiation temperature.", "8.1. GJ 9827 d as a mixed-envelope H 2 O-rich 'steam world'": "The most likely interpretation for the high-metallicity atmosphere with abundant water revealed by all the fits we performed is a 'steam world' mixed atmosphere scenario, with dozens of percent of the atmosphere comprised of high mean-molecular weight species (Figure 4). For the warm equilibrium temperature of GJ 9827 d, we expect its temperature-pressure profile to cross directly from the gas to the supercritical phase of water, without allowing for any water condensation in the atmosphere or at the surface. Since hydrogen and water are highly miscible in the gas and supercritical phases (Soubiran & Militzer 2015; Pierrehumbert 2022; Innes et al. 2023), HMMW volatiles can be expected to be well-mixed throughout the envelope, and the measured upper-atmosphere O/H ratio can be expected to be representative of the deep envelope, with the ratio of H 2 O/CO/CO 2 dictated by the local thermodynamical conditions (Benneke et al. 2024). The measured O/H ratio only provides a lower limit on the oxygen content accumulated by the planet during its formation, as more refractories can be locked in the deep interior (Schlichting & Young 2021; Dorn & Lichtenberg 2021; Luo et al. 2024).", '8.2. nondetection of atmospheric escape in the context of a metal-rich atmosphere': "Beyond the atmospheric composition constraints, our hydrodynamic simulations reveal that a metal-rich atmosphere is favored for GJ 9827 d because of the instability of a hydrogen-dominated atmosphere to atmospheric escape and the nondetection of escaping neutral hydrogen.  \nWhile our hydrodynamic simulations for a hydrogenrich atmosphere scenario predict that the escape would be detectable in Ly α , observations obtained only upper limits (Kasper et al. 2020; Carleo et al. 2021; Krishnamurthy et al. 2023). Furthermore, even at the present-day modeled mass loss rate, the H 2 -He atmosphere would be removed in less than 1 Gyr, while the atmospheric mass loss rates at the beginning of evolution  \nare expected to be 10-100 times higher (Kubyshkina et al. 2024). Therefore, a hydrogen-helium atmosphere is inconsistent with the age of the GJ 9827 system (which is ≥ 1.4 Gyr Passegger et al. 2024).  \nHowever, the mass loss rate is significantly reduced in metal-enriched and steam atmosphere cases. We find that if the atmosphere metallicity is greater than 500 × solar, or for water mass fractions of ≳ 90%, the mass loss is reduced enough to shield the planet from photoevaporative stripping. Moreover, high water mass fractions also lead to an increase by two orders of magnitude in the fraction of hydrogen atoms escaping in the ionized form rather than in the neutral form. This can complicate the detection of hydrogen escape with methods like absorption in the Ly α line, which targets neutral hydrogen.  \nThe escape and atmospheric evolution argument therefore further supports the 'steam world' scenario with a water-rich, high metallicity atmosphere favored by the atmospheric retrieval. In this scenario, and since GJ 9827 d is not a rocky planet from bulk density considerations, the nondetection of an escaping atmosphere has two plausible explanations: the escape rate is too low to be detected (for GJ 9827 d, the expected reduction compared to the H-He dominated atmospheres is up to an order of magnitude), or H and He escape in different states that are not probed by the neutral H (Lyman α ) and metastable He transitions where this search has been conducted to date.", '8.3. Constraints on the interior composition of a water-rich GJ 9827 d': 'We infer a bulk water mass fraction of 32 ± 10% for GJ 9827 d in the case of a water steam atmosphere, which is larger than that of the more massive, larger sub-Neptune TOI-270 d (Benneke et al. 2024). GJ 9827 d would not only have a more metal-rich atmosphere, but it may also have formed from the accretion of a lower relative amount of rocky vs. icy material, as the rock/iron mass fraction inferred for TOI-270 d is ∼ 90%.  \nWe note that our approach could underestimate the total water budget of GJ 9827 d, as significant amounts of water can remain locked in the core or mantle of the planet during mantle crystallization through volatile partitioning(Salvador & Samuel 2023), or could still be dissolved in the magma ocean if at least some of the mantle still is in a molten state (Dorn & Lichtenberg 2021). Another important consideration is the possibility that other volatiles such as CO 2 could be abundant in the planetary atmosphere, making it more compact than this simple model predicts and altering the inferred mass locked in the core and mantle.', '8.4. Volatile enrichment mechanisms': "There are three main pathways to a volatile-rich, and specifically a water-rich atmosphere: accretion of icerich material during or after planet formation, interactions between a magma ocean and the planetary atmosphere, or gradual increase in the atmospheric metallicity from selective loss of lighter species. If GJ 9827 d formed beyond the ice line, it could have accreted large amounts of water ice in solids (L'eger et al. 2004; Luger et al. 2015; Alibert & Benz 2017; Kite & Ford 2018; Venturini et al. 2020; Bitsch et al. 2021). As the planet migrates inward, the photosphere will then become warm enough that a fraction of the accreted water will partition in the vapor form in the atmosphere where it can be probed observationally. Alternatively, magma ocean geochemical interactions with a hydrogen-rich atmosphere can provide an endogenous source for tens of percent of water in the planetary atmosphere, or atmospheric metal mass fractions of ∼ 30 -70% (Schlichting & Young 2021; Kite & Schaefer 2021; Rogers et al. 2024). A third pathway to enhance the metal content of an initially hydrogen-dominated atmosphere which may have been accreted in-situ is the preferential loss of hydrogen over heavier metal species via mass fractionation (Hu et al. 2015; Chen & Rogers 2016; Malsky & Rogers 2020). This process requires a stronger gravitational well for GJ 9827 d than for cooler planets which are more traditionally believed to be in the 'sweet spot' to balance the effect of metals being lost alongside hydrogen. We find that on the ∼ 6 . 2d orbit of GJ 9827 d, its mass of 3.02 M ⊕ is greater than the ≳ 2 . 5 M ⊕ required to undergo mass fractionation in their atmospheres if the dominant mode of atmospheric escape is core-powered mass loss. However, the higher escape fluxes expected from EUV-driven photoevaporation would drive the loss of the bulk of the metal species in the escaping flow (Cherubim et al. 2024). If the temperature is high enough, this can lead to the drag of heavier metals and, therefore, fractionation. Recent research by Louca et al. (in prep.) demonstrated that planets originating as warm super-Neptunes can undergo significant mass loss, shedding both hydrogen and heavier metals. This process results in high-metallicity atmospheres for mature sub-Neptune-like planets. Typically, most of the hydrogen-dominated envelope is lost within the first billion years post-formation, leaving these planets with a considerably lower escape rate in their later stages (see e.g., Lopez et al. 2012; Kubyshkina et al. 2020).", '8.5. A potentially oxidizing atmosphere': 'The high water mixing ratio suggested by our HST /WFC3 + JWST /NIRISS SOSS joint retrieval  \nraises the possibility of a more oxidizing bulk atmosphere composition, which could result from an oxidized mantle composition, the accretion of solid-derived volatiles, or even gradual oxygen buildup over time (Figure 6).  \nOxygen-rich atmospheres have been predicted for temperate water worlds due to the suppression of oxygen sinks by overburden pressure (Krissansen-Totton et al. 2021). They have also been proposed as a possible evolutionary outcome for warm sub-Neptunes interior to the runaway greenhouse limit due to extensive water dissociation and H escape (Harman et al. 2022). Nebular hydrogen in contact with molten silicates will reduce FeO to produce metallic Fe and H 2 O (Kite & Schaefer 2021); the resulting atmosphere may possess a broad range of H 2 :H 2 O inventories depending on nebular endowment, initial mantle FeO, extent of equilibration, and removal of metallic Fe in the core (Schlichting & Young 2022; Kite et al. 2020). In contrast, the accretion of predominantly icy material onto silicates will result in a more oxidized (higher H 2 O:H 2 ratio) bulk atmosphere composition (Kite et al. 2020). For these planets, a plausible amount of H escape could overwhelm mantle oxygen sinks and flip atmospheric composition from net reducing (H 2 O-H 2 dominated) to oxidizing (H 2 O-O 2 dominated) (Harman et al. 2022). The absence of reduced species like CH 4 and NH 3 (Figure A8) and an aerosolfree atmosphere due to the photochemical oxidation of haze-forming species could be consistent with such an atmosphere. The lack of sensitivity to CO/CO 2 over the NIRISS/SOSS wavelength range, however, precludes us from a precise estimation of the total oxygen content in the atmosphere.', '8.6. Evidence for compositional diversity in sub-Neptune atmospheres': "Our characterization of the atmosphere of GJ 9827 d underscores the diversity in the atmospheres of small sub-Neptunes where initial C/N/O/S volatile endowments, (potentially selective) atmospheric escape, and interior-atmosphere interactions over the planet's thermal evolution can all significantly alter the resulting atmosphere. We illustrate this diversity by visually comparing the transmission spectra of TOI-270 d and GJ 9827 d with the larger K2-18 b (Figure 7). While K2-18 b was found to have a hydrogen-dominated, low mean molecular weight atmosphere, with a typical 'gas dwarf' composition (Madhusudhan et al. 2023), the transmission spectra of TOI-270 d and GJ 9827 d both reveal much more metal-enriched atmospheres (Benneke et al. 2024) as reflected in their smaller vertical extent on Figure 7 where they are scaled relative to the number of  \nMantle  \nf  \nO  \nFigure 6. Qualitative description of the major oxygen-, carbon-, nitrogen- and sulfur-bearing species expected in metalrich atmospheres as a function of atmospheric origin and mantle composition (white labels inside of each mantle composition case). Water-dominated atmospheres can result from magma ocean-atmosphere interactions under oxidized mantle conditions (Schlichting & Young 2022; Kite et al. 2020) or from the planet-forming or accreted solids that coagulated beyond the water ice line. Depending on the redox state of the mantle and the H/C ratio of the mantle (Liggins et al. 2021), CH 4 or CO/CO 2 can be the dominant carbon-bearing species, and NH 3 or N 2 can be the dominant nitrogen-bearing species. Further O 2 buildup can occur over time as H is gradually lost to space (Krissansen-Totton et al. 2021; Harman et al. 2022).  \n<!-- image -->  \nplanetary scale heights spanned by the spectral features (calculated for a typical H 2 /He dominated atmosphere).  \nThe trend of decreasing, then increasing spectral feature strength for planets with warm-temperate equilibrium temperatures, with a 'minimum' found near ∼ 600 K (Brande et al. 2023) is also reflected in the gradual decrease in feature strength from the colder K218 b to the warmer TOI-270 d and GJ 9827 d. However, while initially attributed to the effect of verticallyextended clouds, or small-particle hazes as a function of equilibrium temperature (Brande et al. 2023), clouds were not significantly detected in the atmospheres of TOI-270 d or GJ 9827 d, with the variety of spectral feature strengths seemingly arising from differences in atmospheric mean molecular weight, rather than cloudiness trends. This motivates further exploration of the diverse atmospheres of small sub-Neptunes to reveal potential trends that, beyond equilibrium temperature, could track with planet mass, stellar type, or planetary surface gravity.", '8.7. Observational discriminants of atmospheric origins': 'Further observational probes of the atmosphere of GJ 9827 d have the potential to shed more light on its HMMW volatile makeup and on the origins of its volatile-rich atmosphere. In particular, upcoming transit observations with NIRSpec/G395H as part of GO  \n4098 will open the possibility of detecting carbon species such as CO and CO 2 , or even sulfur-bearing species (e.g. SO 2 ). The constraining power of NIRISS/SOSS for CO or CO 2 abundances is very limited (see Figures A9, A10, Table 2), while both these molecules have larger cross sections and major absorption features in the NIRSpec/G395H bandpass.  \nA detection of CO 2 as the dominant carbon-bearing species would confirm our tentative conclusion of the oxidizing state of the atmosphere and suggest more C than H in the initial silicate endowment of the mantle. However, more detailed modeling would be required to tell apart the origin of the oxygen, either from outgassing from a highly oxidized mantle or direct volatile enrichment from the protoplanetary disk (Figure 6). Abundant SO 2 is a strong signature of oxidized mantle outgassing, but can also be produced by upper atmosphere photochemistry of H 2 S and H 2 O (Tsai et al. 2023). Alternatively, if CO dominates, the favored interpretation for the volatile enrichment would be mantle outgassing rather than the accretion of solid-derived volatiles as CO 2 would quickly dominate in the most oxidizing conditions. Finally, the detection of abundant methane or ammonia, which are only weakly constrained by our NIRISS/SOSS spectrum, would suggest a carbon-poor silicate mantle (Figure 6).  \nFurther, the atmospheric CO 2 /CH 4 ratio was recently proposed as an observational window into the O/H or  \nFigure 7. Transmission spectra of K2-18 b (Madhusudhan et al. 2023), TOI-270 d (Benneke et al. 2024), and GJ 9827 d, in units of scale heights for a H 2 /He-dominated atmosphere, offset. For clarity, the spectra of TOI-270 d and GJ 9827 d are binned to a resolving power of 50 (gray points). The bestfit models from the retrievals run on the spectra of GJ 9827 d and TOI-270 d are shown on top of the observations (colored lines) and for K2-18 b we show a representative SCARLET model for the retrieved abundances quoted in Madhusudhan et al. (2023). The background color and line colors encode the variation in equilibrium temperature between the three sub-Neptunes, but they also have different sizes and densities (see annotated temperatures for a Bond albedo of 0.3, masses, and radii).  \n<!-- image -->  \nH 2 O/H 2 ratio in the deep interior of warm-temperate planets (Yang & Hu 2024). If this finding holds for the ∼ 620 K GJ 9827 d, a detection of carbon-bearing species could help us differentiate between water derived from the accretion of ice-rich solids, or from volatile-poor formation followed by interior atmosphere interactions, although in both cases volatile partitioning between the core, mantle, and atmosphere will alter the atmospheric volatile content. Observations at longer wavelengths could potentially reveal oxygen buildup via the presence of O 2 (via the collisionally induced absorption band at ∼ 6 µ m; Fauchez et al. 2019) or O 3 at 9.6 µ m, but further theoretical work is required to estimate the expected abundance of these species in the atmospheres of small sub-Neptunes such as GJ 9827 d and assess their detectability.', '9. SUMMARY AND CONCLUSIONS': "We obtained two new transit observations of GJ 9827 d with JWST/NIRISS SOSS as part of the JWST subNeptune survey (GO 4098). We detect seven candidate flare events in the data, three of which are confirmed as bona fide flares by their time evolution and measured spectral energy distribution. The spectra obtained from both visits are best explained by water features with small ∼ 20 -40ppm amplitudes. Our water detection, in agreement with the HST/WFC3 result (Roy et al. 2023), is robust even when the TLS effect is considered.  \nThe addition of the JWST/NIRISS SOSS observations enables us to break the degeneracy between lowmetallicity cloudy and high mean molecular weight atmospheres, and we detect a high mean molecular weight, metal-enriched and potentially water-rich 'steam world' atmosphere which from irradiation considerations would be in a well-mixed vapor/supercritical state. The metalenriched atmosphere scenario is also supported by our and previous nondetections of escaping neutral hydrogen and metastable helium from the atmosphere of GJ 9827 d as escape rates are lower in higher mean molecular weight atmospheres and hydrogen ionization fraction in the escaping flow is larger, eluding observational probes. We estimate, for a pure-H 2 O envelope on top of a solid rock/iron core, that the hydrosphere makes up between 20 and 40% of the total planet mass.  \nThe atmospheric metal enrichment could originate from large initial volatile inventories from early accretion, preferential loss of H through mass fractionation, or enrichment through magma ocean-atmosphere interactions over geological timescales. Future observations probing potential absorption from carbon- and sulfurbearing species could shed more light on which volatile enrichment mechanism(s) shaped the atmosphere of GJ 9827 d.  \nWe thank the anonymous referee for providing comments that improved our manuscript. This work is based on observations with the NASA/ESA/CSA James Webb Space Telescope, obtained at the Space Telescope Science Institute (STScI) operated by AURA, Inc. All of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute. The data used in this paper can be found in MAST at the following DOIs: 10.17909/sx19-8852 ( JWST ) and 10.17909/dvqh-2r48 ( HST ). The authors thank O. Lim for providing a NIRISS/SOSS throughput file and A. Feinstein for sharing her perspective on stellar flares and their impact on this dataset. C.P.-G. acknowledges support from the NSERC Vanier scholarship and the Trottier Family Foundation. This work was funded by the Institut Trottier de Recherche sur les Exoplan'etes (iREx). M.R. acknowledges funding from NSERC, Fonds de recherche du Qu'ebec - Nature et technologies (FRQNT), and iREx. R.A. acknowledges the SNSF support under the Post-Doc Mobility grant P500PT 222212 and the support of the Institut Trottier de Recherche sur les Exoplan'etes (iREx). M.F.T. acknowledges financial support from the Clarendon Fund and the FRQNT.  \nTable A1. Fitted planetary parameters to both JWST /NIRISS SOSS Order 1+2 transits of GJ 9827d extracted using supreme-SPOON and fitted using ExoTEP . We report the planet-to-star radius ratios separately for the white-light-curve fits to each order.  \nFigure A1. Comparison of the fitted transit time variations for visits 1 and 2 with NIRISS/SOSS (black points) to the previous measurements from Kepler (blue), HST (purple), and Spitzer (pink). We show the transit-timing variations (TTVs) as a function of epoch, after fitting a new ephemeris to all the transit times combined. The inset zooms in on the two NIRISS/SOSS transits presented in this work. The JWST transits nearly double the baseline, occurring over 1000 days after the last HST transit, improving the precision of the orbital period measurement by a factor of five.  \n<!-- image -->", 'A. COMPARISON OF FITTED PARAMETERS WITH LITERATURE VALUES': 'Our two measured transit times enable us to refine the ephemeris of GJ 9827 d (see Table 1). Not only is the precision of our transit times measurement higher compared to previous sources, but the two NIRISS/SOSS transits also occurred over 1000 days after the previously-measured HST transit times (Figure A1, Roy et al. 2023).  \nWe do not detect significant transit-timing variations from the two visits, separated by only one planetary orbit. This is in contrast to the 10 measured HST transit times, which exhibit TTVs on the order of 5 to 10 minutes (Figure A1). We note, however, that the HST transit times might be biased by the assumed b and a/R ⋆ fixed in the HST lightcurve fits to literature values (Niraula et al. 2017) which differ significantly from our retrieved impact parameter and semi-major axis (Table A1). Our fitted orbital parameters and the derived planetary properties agree with literature values (Kosiarek et al. 2021; Passegger et al. 2024). Our derived planetary radius has larger uncertainties than the Spitzer value even if our R p /R ⋆ is more precise, with the error budget dominated by the larger error bars on the stellar radius (0 . 58 ± 0 . 03 R ⊙ ) from Passegger et al. (2024) compared to the isochrone determination from Kosiarek et al. (2021) (0 . 579 ± 0 . 018 R ⊙ ). The more conservative error bars from Passegger et al. (2024) are derived using a joint GP modeling of the transits and radial velocity data with a broad prior informed not only by the isochrone fit to the ESPRESSO-derived stellar parameters but also other determinations using alternative age determination methods.', 'B. CONSISTENCY OF SPECTRUM ACROSS BINNING SCHEMES, DATA REDUCTIONS, AND INSTRUMENTS': 'We compare for each visit the spectra extracted 1- and 2-pixel resolution with the R=100 spectrum for the supreme-SPOON reduction and find that they agree well within 1 σ over the wavelength range they share in com-  \nFigure A2. Comparison between the spectra that were extracted at high (1-pixel in blue and 2-pixel in orange) resolution, and the R=100 (black) spectrum, for Visit 1 (top panel) and Visit 2 (bottom panel) using the supreme-SPOON version of the data reduction. The three versions of the transmission spectrum for each visit are consistent to well within their quoted uncertainties in each spectral bin.  \n<!-- image -->  \non (Figure A2). We also compare the spectra obtained from the supreme-SPOON and NAMELESS reductions for each visit, and find a similarly satisfactory agreement (Figure A3). Finally, we compare the JWST spectrum to the HST /WFC3 transmission spectrum of GJ 9827d and find that they match well over the bandpass of HST /WFC3, except for a small vertical shift (Figure A4). The transmission spectrum is provided in Table A2 and in a machinereadable format along with the published article.', 'C. ANALYSIS OF CANDIDATE STELLAR FLARES': "The multiple long and short timescale events in the observed light-curves can naturally be explained by stellar activity, which principally consists of stellar flares or spot-crossing events. Both events will display different temporal and spectral variations and should be distinguishable in our data. On one hand, the red optical/NIR spectrum of flares evolves on timescales of ∼ 10 s. It is composed of a rapid period of prompt emission with an effective temperature ( T eff ) of ∼ 9000 K followed by a more gradual period of emission with a T eff of ∼ 5000 K (Kowalski et al. 2016). As a result, the observed T eff of the post-subtraction stellar spectrum rapidly increases and then slowly decays during a flare and follows a fast-rise, exponential decay (FRED) temporal profile. Conversely, variations in the residual spectrum resulting from a spot-crossing event occur more slowly and with a symmetrical temporal profile well-described by a Gaussian (e.g., B'eky et al. 2014; Schutte et al. 2023; Biagiotti et al. 2024). Furthermore, a flare spectrum is brightest at near-UV and blue optical wavelengths (0.2-0.5 µ m), placing the NIRISS/SOSS wavelength range in the Rayleigh-Jeans tail during the flare peak. Meanwhile, starspot spectra are generally brighter at IR than optical wavelengths (Waalkes et al. 2024; Biagiotti et al. 2024). Although it only enables us to sparsely sample the potentially rapidly-evolving brightness profiles of flares, the 16.5-s cadence of the NIRISS/SOSS observations is deemed sufficient to determine whether the T eff evolution favors the interpretation of a mid-transit event as either a spot crossing event or flare through the maximum T eff value and the presence of (a)symmetry in the T eff evolution.  \nWe flux-calibrate the extracted spectra from the supreme-SPOON pipeline against a K7V PHOENIX spectrum obtained from the MUSCLES Extension archive (Behr et al. 2023), and follow the procedure applied to NIRISS/SOSS observations of TRAPPIST-1 (Howard et al. 2023). The K7V PHOENIX spectrum is computed from the Lyon  \nFigure A3. Comparison between the spectra obtained from the supreme-SPOON (blue) and NAMELESS (orange) data reductions. The top panel shows the R=100 spectrum for visit 1, and the bottom panel shows the same thing for visit 2. The spectra agree to better than 1 sigma in almost all spectral bins.  \n<!-- image -->  \nBT-Settl CIFIST 2011 2015 grid (Allard 2016) for a star of T eff =4124 K, d ∗ =86.7 ± 0.3 pc, M ∗ =0.72 ± 0.03 M ⊙ , and R ∗ =0.67 ± 0.01 R ⊙ , and scaled to match the distance to GJ 9827. Comparison with a blackbody spectrum at the effective temperature of GJ 9827 reveals that the flux in the atmosphere model exceeds the blackbody flux by a factor of 1.348 in the wavelength range of 0.6-1.6 µ m where the majority of flare emission occurs in NIRISS observations from the shape of the flare blackbody spectrum. As a result, we further scale the PHOENIX spectrum by a factor of 0.742 to maximize calibration accuracy within this region. The color correction for each visit is obtained by taking the ratio between the PHOENIX model spectrum and an average spectrum of all quiescent integrations during the visit. For the construction of this quiescent spectrum, we select integrations that are both out of transit and in time intervals where the stellar activity (assessed from the NIRISS/SOSS light-curve integrated in the TESS bandpass) is low. Before computing the ratio of the observed and model spectra, we bin the spectrum in resolution (200 columns per wavelength bin for order 1 and 160 for order 2) in order to remove localized spectral features (Howard et al. 2023). This ratio is then used to flux-calibrate each of the extracted spectra in the time series and to obtain a flux-calibrated quiescent spectrum.  \nWe construct flare-only spectra by calculating the difference between each flux-calibrated spectrum in the time series (normalized by the transit light-curve), and the quiescent spectrum. To obtain the batman model shown in Figure A5, we choose to fix the transit depth in the model fitted to the visit 2 light-curve to that obtained from a fit to the visit 1 transit because of the long decay tail of the largest F3 flare. This causes a difference in apparent 'fit quality' with the white light curve fit used to obtain the transmission spectrum (Figure 2) but ensures that the flare analysis is not biased by to degeneracy between the flare tail and transit depth. We carefully select the lowest-activity epochs when computing the median out of transit flux. The measured flare properties are more sensitive to the inclusion of small flux excursions in the quiescent baseline than are the properties of a time-averaged transmission spectrum, requiring a separate batman model subtraction.  \nFigure A4. Comparison of the spectrum of GJ 9827 d over both JWST/NIRISS SOSS transit epochs and with the HST/WFC3 observations. Top panel: Transmission of GJ 9827 d as observed during the first (blue) and second (orange) visit, with the spectrum of visit 1 shifted down by about 40 ppm to match the median value of the visit 2 spectrum. The combined spectrum from both visits is shown in black. Bottom panel: Combined JWST/NIRISS SOSS transmission spectrum (gray) along with binned version (black) and the spectrum as observed from 10 HST/WFC3 transits (Roy et al. 2023), shifted up by 10 ppm.  \n<!-- image -->  \nWe fit the T eff of the flare-only spectrum for each integration using a single-temperature Planck function model (Howard et al. 2023), allowing us to obtain a flare effective temperature time series for each visit (Figure A5. We also use the convolution of each NIRISS/SOSS spectrum with the transmission curve of TESS in the 0.6-1.0 µ m range, to construct a light-curve that can be compared with the light-curves of thousands of K and M dwarf flares observed at a comparable 20-s cadence by the TESS mission (e.g. Gilbert et al. 2022; Howard & MacGregor 2022; Paudel et al. 2024).  \nWe flag as flares events in the post-processed TESS band light-curve where at least two consecutive integrations exceed the measured visit-specific scatter by ≥ 3 σ . Flare energies are calculated using the quiescent energy of GJ 9827 in the TESS bandpass ( Q 0 =1.000 ± 0.002 × 10 32 erg s -1 ) and the flare equivalent duration (ED; Howard et al. 2023). We detect seven small candidate flare events with TESS band energies ranging from 7.40 ± 0.60 × 10 30 to 3.85 ± 0.12 × 10 31 erg across both visits (labeled F1-F7; see Figure A5). The presence of seven flares with E TESS ≥ 7.4 × 10 30 erg during a total stare time of 0.289 d qualifies GJ 9827 as an active K-dwarf, whose flare rate is typical of other active field stars of spectral class K7-M0 (Howard et al. 2019).  \nWe measure the T eff values from the flare-only spectra of each event, averaging over all integrations during the peak times of the flare to increase the signal to noise. The S/N of the T eff time series of the F1, F3, and F4 events are sufficient to observe the rapid rise and slow decay profile typical of stellar flares, but not starspots. For example, the rise time of the F1, F3, and F4 T eff time series lasts 2.17, 0.29, and 0.27 min, respectively. The rise times of spot crossing events during transits of HAT-P-18 b and HATS-2 b are 6.9, 11.2, and 13.6 min, respectively (Fournier-Tondreau et al. 2023; Biagiotti et al. 2024). The T eff value obtained from the flare-only spectrum averaged over the peak integrations for each of these same events is 4470 ± 320, 4580 ± 230, and 3980 ± 360 K, respectively (Figure A6).  \nThe peak continuum temperatures are similar to those of four flares from the M8 dwarf TRAPPIST-1 observed with NIRISS and NIRSpec. Although for GJ 9827 these temperatures are close to the effective temperature of the star, the flux amplitude of the flare-only spectrum used to measure the T eff values is ≥ 5 times larger than the uncertainty in the subtraction of the quiescent stellar spectrum from 0.6-1.7 µ m, and ≥ 50 larger than the uncertainty at 0.9 µ m. The TRAPPIST-1 flares also display a FRED time series in T eff that peaks at values of up to 5300 K (Howard et al. 2023). We fit an empirical stellar flare template model (Tovar Mendoza et al. 2022) and a Gaussian model to the T eff time series of flares F1, F3, and F4 and measure the residual sum of squares (RSS) for each model to quantify the degree of asymmetry. The flare model is preferred in each case by a factor of 3.1, > 5, and 4.1 for F1, F3, and F4, respectively. We conclude that the combination of the extremely rapid and asymmetrical continuum temperature evolution and the 4000-5000 K peak T eff values of these events are best explained by the flare interpretation.  \nWhile for TRAPPIST-1, the H α line exhibits variations that track the time evolution of the flare, we do not find clear flare evidence in H α for GJ 9827. Although the TESS-band flare energies we obtain for GJ 9827 are similar to those of TRAPPIST-1, a quiescent K7V star is three orders of magnitude brighter at 0.656 µ m than an M7.5V star, such that a flare from GJ 9827 must be 10 3 times more luminous in the H α line to produce the same H α light-curve as a flare from TRAPPIST-1. The TRAPPIST-1 flares emitted only 10% as much energy in the H α line as in the continuum. If this energy distribution holds for late K dwarfs, the presence of continuum enhancements is, therefore, more indicative of flaring than are H α enhancements, which could explain our nondetection.  \nWe use optical-NUV and NUV-FUV flare energy scaling relations to translate our measured flare energies to FUV (0.09-0.17 µ m) and NUV (0.17-0.32 µ m) energies (Table A3). Even small flares observed at NIR wavelengths can produce significant enhancements in the stellar emission at UV wavelengths due to the peak of the flare SED from 0.2-0.4 µ m (Jackman et al. 2023; Paudel et al. 2024). We employ the E TESS , FWHM -E NUV , tot relation obtained from simultaneous 20-s cadence TESS and Swift observations to estimate the total NUV energy emitted during the F1-F7 flares (Howard et al., in preparation). This relation leverages the time resolution of the 20-s cadence TESS data to measure the energy emitted in the TESS band with a two-fold increase in precision. We propagate the mean and 1 σ uncertainty in the estimated NUV energies into the FUV using the ratio of flare emission at FUV to NUV wavelengths. This ratio, R FUV / NUV = 0 . 49 ± 0 . 02, was recently measured from GALEX observations of a sample of 182 flares from 158 field stars which provide a strong constraint on the FUV emission of any NUV flare (Berger et al. 2023).  \nThe combined FUV and NUV energies of F1-F7 make up 32% and 3% of the integrated quiescent stellar FUV and NUV stellar emission throughout all transit observations, respectively. The flares can contribute larger fractions of the total UV irradiation over shorter time intervals. For example, the F3 event contributes 82% and 8% of the integrated quiescent emission at FUV and NUV wavelengths respectively during the first hour of the transit observation during visit 2. Variable stellar emission at FUV and NUV wavelengths can impact the abundances of volatile species in the upper atmosphere through photochemical processes, and the frequent flares we observe suggest that photodissociation might be even more prevalent (especially during the second NIRISS/SOSS transit) than estimated from observations of the quiet state of GJ 9827.", 'D. MODELING OF THE TRANSIT LIGHT SOURCE EFFECT': 'If the stellar surface is heterogeneous, in the presence of dark spots or bright faculae, the average stellar surface spectrum may not be a good representation of the spectrum of the star along the transit chord (see e.g. Sing et al. 2011; Berta et al. 2011; Rackham et al. 2018). In particular, for the case of GJ 9827 d, unocculted stellar spots can mimic water features in the transmission spectra of planets around cool stars (e.g., Barclay et al. 2021; Brande et al. 2022).  \nIn that case, the transmission spectrum of the planet is multiplied by a wavelength-dependent factor ϵ λ, het (Rackham et al. 2018; Pinhas et al. 2018):  \nϵ λ, het = 1 1 -N het ∑ i =1 f het , i ( 1 -I λ, het , i ( T het , i ) I λ, phot ( T phot ) ) (D1)  \nwhere f het , i describes the fraction of the stellar disk covered by heterogeneity i , while I λ, het , i , and I λ, phot correspond respectively to the specific intensity of the stellar heterogeneity i with temperature T het , i and of the stellar photosphere with temperature T phot .  \nIn our modeling, we assume one or two populations of stellar heterogeneities:  \nFlux change [ppt]  \nFigure A5. Top panels : WL light-curves showing the full baseline versus the transit model. The 1 σ flux range during times of quiescence is highlighted. We opt to use the transit depth from visit 1 during visit 2 due to degeneracy between the depth and the flare decay tail of F3. Bottom left panels : The WL transit light-curve, TESS band time series, and T eff time series for the flare-only continuum spectra are shown for visit 1. Times during flares are highlighted, where flares are defined as events with TESS band flux increases with two or more points exceeding the local noise by 3 σ in the binned light-curve. The 16.5-s cadence data are shown in gray, while binned time series are shown in blue and black for the TESS and T eff time series, respectively. A fitted linear trend in the out-of-transit light-curve has been removed from the WL light-curve to better show flare events. The transit model shown in the top panel has been subtracted from the TESS time series. Gaps in the T eff time series are due to low S/N or nondetections outside of the flare peak. Bottom right panels: Same as the bottom left panels, but for visit 2.  \n<!-- image -->  \nFigure A6. Time-averaged spectra of all integrations during the flare peak times of the F1-F4 candidates. Each spectrum is obtained by subtracting the out-of-flare stellar spectrum to create a flare-only spectrum. The flare-only spectrum is binned in wavelength and fit with a single-temperature blackbody model to obtain the hottest T eff reached during each event.  \n<!-- image -->  \n- 1. A single population of stellar spots colder than the stellar photosphere: N het = 1 in Eq. D1 and T spot < T phot ;\n- 2. Two populations of heterogeneities ( N het = 2 in Eq. D1) corresponding to stellar spots ( T spot < T phot ) and faculae ( T fac > T phot ), respectively.  \nIn our modeling, the spot temperature is parameterized as T spot = T phot +∆ T spot , with ∆ T spot < 0, and ∆ T fac > 0 similarly describes the temperature of the faculae component in the two-heterogeneities case.  \nWe display the retrieved distribution of stellar contamination-only models that provide the best match to the NIRISS/SOSS + HST/WFC3 spectrum of GJ 9827 d in the spots-only case, and find that the observed water features are inconsistent with stellar contamination alone (Figure A7).', 'E. CHAIN ATMOSPHERIC ESCAPE MODELING SETUP': 'The CHAIN model combines a 1D hydrodynamic upper atmosphere model based on that of Kubyshkina et al. (2018) solving the hydrodynamic equations with the non-local thermodynamical equilibrium (NLTE) photoionization and radiative transfer code Cloudy (Ferland et al. 2017) providing detailed (photo)chemistry calculations, hence, realistic heating and cooling rates. The Cloudy code can account for the 30 lightest chemical elements from hydrogen to zinc and supports various compositions (namely, abundances of individual elements (He-Zn) can be set to zero or scaled relative to hydrogen); the code can also account for a range of common molecular species, which are typically present below the 1 µ bar pressure level and destroyed in photodissociation reactions at higher altitudes (e.g., Egger et al. 2024). The hydrodynamic part of the CHAIN code assumes that the atmosphere is well-mixed within the simulation domain (for close-in planets typically set between 1-100 mbar level and 1-2 Roche radii of the planet) and treats it as a uniform outflow with the mean molecular weight adjusted to reproduce the employed composition. To verify this assumption, the code controls that the exobase position is above the sonic point throughout the simulation.  \nThe code accounts for the full stellar energy distribution between infrared and X-ray. As a proxy for GJ 9827, we employ the shape of the spectrum of WASP-43 provided by the Measurements of the Ultraviolet Spectral Characteristics of Low-mass Exoplanetary Systems (MUSCLES) Treasury Survey (version 23, France et al. 2016; Youngblood et al. 2016; Loyd et al. 2016. WASP-43 b has a mass, temperature, and luminosity similar to those of GJ 9827. We further scale the visible-infrared part of the spectrum to reproduce the bolometric luminosity of GJ 9827 and the X-ray plus  \nTable A2. Truncated NIRISS/SOSS spectrum of GJ 9827 d, using the supreme-SPOON reduction, fitted at R=100. A full version of this spectrum is provided as a machine-readable table in the electronic journal version.  \nTable A3. Properties of flare candidates from GJ 9827 during visits 1 and 2. Columns are flare ID, flare start and stop time, peak significance of detection in the binned TESS band light-curve, TESS band flare energy, estimated NUV flare energy, estimated FUV energy, and T eff of the flare-only spectrum at peak. The NUV and FUV energies are derived from optical-NUV scaling relations based on simultaneous Swift and TESS flare observations in Paudel et al. (2024) and Howard et al. (in prep.), and GALEX flare observations in Berger et al. (2023).  \nextreme ultraviolet part (XUV) and anchor our results in the previous literature estimate (Carleo et al. 2021), who calculate the XUV flux received by GJ 9827 d of ∼ 4300 erg s -1 cm -2 using the method outlined in Linsky et al. (2014). This value could potentially be larger during flare events (Section C). We further verified that the shape of the spectrum does not have a crucial effect on our results by re-running the simulations employing Sun-like and GJ 436 spectra; we found that it affects to some extent the abundances of ion species at low altitudes but not the escape rates. The overall effects are similar to those described in Kubyshkina et al. (2024). We note that CHAIN does not employ a multi-fluid approach, and we rely on the estimations of the maximum mass of species dragged by the hydrodynamic  \nFigure A7. Result from the stctm fit (spots only) to the NIRISS/SOSS + HST/WFC3 spectrum of GJ 9827 d. The combined JWST+HST spectrum is binned to R=50 for visual purposes (black points). The 1, 2, and 3 σ constraints on the retrieved stellar contamination-only spectra are shown (semi-transparent regions) as well as the median (coral line) and best-fit spectrum (black line). Stellar contamination without a planetary atmosphere signal cannot reproduce the water features observed in the transmission spectrum of GJ 9827 d (see Figure 3 for an atmosphere+stellar contamination fit).  \n<!-- image -->  \noutflow to (Hunten et al. 1987), which may suffer from some inconsistencies (Linssen et al. 2024). Therefore, it is not clear at this stage whether the oxygen ions present at high altitudes would escape along with the hydrogen.', 'F. MODEL COMPARISON AND POSTERIOR DISTRIBUTIONS FROM ATMOSPHERE RETRIEVALS': 'We use Bayesian model comparison to test the significance of the detection of stellar heterogeneities and H 2 O in the JWST NIRISS/SOSS + HST/WFC3 spectrum of GJ 9827 d (Table A4). We also display the posterior distributions on H 2 O, CH 4 , NH 3 as well as stellar contamination parameters for various retrieval scenarios, with or without stellar contamination and for different assumptions about the background gas in the prior, to illustrate the robustness of our results across retrieval assumptions and frameworks (Figure A8). The posterior distributions on all the fitted molecules are also provided for the SCARLET and POSEIDON retrievals (Figures A9 and A10).', 'G. INTERNAL STRUCTURE MODELING FOR GJ 9827D': 'We show the posterior distribution on all the parameters fitted in the Bayesian framework smint over the grid of water-rich models, which provide us a constraint on the water mass fraction of GJ 9827 d (Figure A11).', 'REFERENCES': "Aguichine, A., Mousis, O., Deleuil, M., & Marcq, E. 2021, The Astrophysical Journal, 914, 84,  \ndoi: 10.3847/1538-4357/abfa99 Albert, L., Lafreni'ere, D., Ren'e, D., et al. 2023, Publications of the Astronomical Society of the Pacific, 135, 075001, doi: 10.1088/1538-3873/acd7a3 Alibert, Y., & Benz, W. 2017, Astronomy & Astrophysics, 598, L5, doi: 10.1051/0004-6361/201629671  \nAllard, F. 2016, in SF2A-2016: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, ed. C. Reyl'e, J. Richard, L. Cambr'esy,  \nM. Deleuil, E. P'econtal, L. Tresse, & I. Vauglin, 223-227  \nAzzam, A. A. A., Yurchenko, S. N., Tennyson, J., & Naumenko, O. V. 2016, Mon. Not. R. Astron. Soc., 460, 4063, doi: 10.1093/mnras/stw1133  \nBarclay, T., Kostov, V. B., Col´on, K. D., et al. 2021, The Astronomical Journal, 162, 300,  \ndoi: 10.3847/1538-3881/ac2824  \nBehr, P. R., France, K., Brown, A., et al. 2023, AJ, 166, 35, doi: 10.3847/1538-3881/acdb70  \nB´eky, B., Kipping, D. M., & Holman, M. J. 2014, MNRAS, 442, 3686, doi: 10.1093/mnras/stu1061  \nBenneke, B. 2015, arXiv e-prints, 1504, arXiv:1504.07655.  \nhttp://adsabs.harvard.edu/abs/2015arXiv150407655B  \nTable A4. Retrieval results comparison metrics. We list model comparison statistics: the log Bayesian evidence ln Z , and the ∆ln Z for each nested model assumption compared to the baseline model without this assumption. For both the SCARLET and POSEIDON retrievals with H 2 /He as the filler gas, the other molecules included were H 2 O (except when explicitly removed), N 2 , CH 4 , NH 3 , CO, CO 2 , HCN, and H 2 S. The results with the centered-log-ratios prior were obtained using POSEIDON. The 'sigma' confidence levels quoted in the right-most column are calculated from the log-evidence values using the standard Bayesian comparison scale (Trotta 2008; Benneke & Seager 2013). The chemically-consistent results are quoted for SCARLET retrievals without stellar contamination.  \nBenneke, B., & Seager, S. 2012, The Astrophysical Journal,  \n753, 100, doi: 10.1088/0004-637X/753/2/100 -. 2013, The Astrophysical Journal, 778, 153, doi: 10.1088/0004-637X/778/2/153 Benneke, B., Werner, M., Petigura, E., et al. 2017, The Astrophysical Journal, 834, 187, doi: 10.3847/1538-4357/834/2/187 Benneke, B., Knutson, H. A., Lothringer, J., et al. 2019a, Nature Astronomy, 3, 813, doi: 10.1038/s41550-019-0800-5 Benneke, B., Wong, I., Piaulet, C., et al. 2019b, The Astrophysical Journal Letters, 887, L14, doi: 10.3847/2041-8213/ab59dc Benneke, B., Roy, P.-A., Coulombe, L.-P., et al. 2024, JWST Reveals CH $ 4 $ , CO $ 2 $ , and H $ 2 $ O in a Metal-rich Miscible Atmosphere on a Two-Earth-Radius Exoplanet, arXiv. http://arxiv.org/abs/2403.03325 Berger, V. L., Hinkle, J. T., Tucker, M. A., et al. 2023, arXiv e-prints, arXiv:2312.12511, doi: 10.48550/arXiv.2312.12511  \nBerta, Z. K., Charbonneau, D., Bean, J., et al. 2011, The Astrophysical Journal, 736, 12,  \ndoi: 10.1088/0004-637X/736/1/12 Biagiotti, F., Mancini, L., Southworth, J., et al. 2024, A&A, 685, A131, doi: 10.1051/0004-6361/202347872 Bitsch, B., Raymond, S. N., Buchhave, L. A., et al. 2021, arXiv e-prints, 2104, arXiv:2104.11631. http://adsabs.harvard.edu/abs/2021arXiv210411631B Bonomo, A. S., Dumusque, X., Massa, A., et al. 2023, Astronomy and Astrophysics, 677, A33, doi: 10.1051/0004-6361/202346211 Brande, J., Crossfield, I. J. M., Kreidberg, L., et al. 2022, The Astronomical Journal, 164, 197, doi: 10.3847/1538-3881/ac8b7e -. 2023, Clouds and Clarity: Revisiting Atmospheric Feature Trends in Neptune-size Exoplanets, arXiv.  \nhttp://arxiv.org/abs/2310.07714  \nBuchner, J., Georgakakis, A., Nandra, K., et al. 2014, Astronomy and Astrophysics, 564, A125, doi: 10.1051/0004-6361/201322971  \nFigure A8. Posterior distributions for important atmospheric volatile species (top, blue) and the stellar heterogeneity component (bottom, orange), for the POSEIDON retrieval (top rows) using a centered-log-ratios (CLR) agnostic prior on the background gas, and for the SCARLET retrieval (bottom rows) where H 2 /He is the filler (background) gas. Different colors correspond to different retrievals where spots or faculae are included or excluded (see legend). Our constraints on the HMMW atmospheric volatile species are broadly unaffected by our treatment of stellar contamination, and the retrieval only provides a lower limit on the H 2 O abundance. For the SCARLET retrieval, the distributions on T spot and T fac are derived from the posterior distributions on ∆ T spot , ∆ T fac and T phot , star . The POSEIDON prior on the spot temperature extends all the way to 2300 K, while SCARLET assumes that spots are at most 800 K colder than the photosphere.  \n<!-- image -->  \nBurn, R., Mordasini, C., Mishra, L., et al. 2024, A radius valley between migrated steam worlds and evaporated rocky cores, arXiv. http://arxiv.org/abs/2401.04380  \nBushouse, H., Eisenhamer, J., Dencheva, N., et al. 2023,  \nJWST Calibration Pipeline, 1.12.3, Zenodo,  \ndoi: 10.5281/zenodo.8404029  \nCarleo, I., Youngblood, A., Redfield, S., et al. 2021, The Astronomical Journal, 161, 136,  \ndoi: 10.3847/1538-3881/abdb2f  \nChen, H., & Rogers, L. A. 2016, The Astrophysical Journal, 831, 180, doi: 10.3847/0004-637X/831/2/180  \nCherubim, C., Wordsworth, R., Hu, R., & Shkolnik, E. 2024, Strong fractionation of deuterium and helium in sub-Neptune atmospheres along the radius valley, arXiv. http://arxiv.org/abs/2402.10690  \nColes, P. A., , Yurchenko, S. N., & Tennyson, J. 2019, Mon. Not. R. Astron. Soc., 490, 4638, doi: 10.1093/mnras/stz2778  \nCoulombe, L.-P., Benneke, B., Challener, R., et al. 2023, Nature, 620, 292, doi: 10.1038/s41586-023-06230-1  \nDarveau-Bernier, A., Albert, L., Talens, G. J., et al. 2022,  \nPublications of the Astronomical Society of the Pacific,  \n134, 094502, doi: 10.1088/1538-3873/ac8a77  \nDorn, C., & Lichtenberg, T. 2021, arXiv:2110.15069 [astro-ph]. http://arxiv.org/abs/2110.15069  \nFigure A9. Posterior distributions on the molecular abundances for all the fitted molecules in the SCARLET retrieval with H 2 /He used as a filler gas and without stellar contamination. The 1, 2, and 3 σ lower (upper) limits for H 2 O (other molecules) are indicated by the color shadings.  \n<!-- image -->  \nFigure A10. Posterior distributions on the molecular abundances for all the fitted species in the POSEIDON retrieval with centered-log-ratios priors and without stellar contamination. The 1, 2, and 3 σ lower (upper) limits for H 2 O (other molecules) are indicated by the color shadings.  \n<!-- image -->  \nFigure A11. Composition of GJ 9827 d for a 100% water envelope. Joint and marginalized posterior distributions of the planet structure fit for GJ 9827 d assuming a vapor/supercritical hydrosphere on top of an interior composed of various mass fractions of rock and iron. The 1, 2, and 3 σ probability contours are shown.  \n<!-- image -->  \n- Salvador, A., & Samuel, H. 2023, Icarus, 390, 115265, doi: 10.1016/j.icarus.2022.115265\n- Schlichting, H. E., & Young, E. D. 2021, Chemical equilibrium between Cores, Mantles, and Atmospheres of Super-Earths and Sub-Neptunes, and Implications for their Compositions, Interiors and Evolution, Tech. rep. https://ui.adsabs.harvard.edu/abs/2021arXiv210710405S\n- -. 2022, The Planetary Science Journal, 3, 127, doi: 10.3847/PSJ/ac68e6\n- Schutte, M. C., Hebb, L., Wisniewski, J. P., et al. 2023, AJ, 166, 92, doi: 10.3847/1538-3881/ace59c\n- Sing, D. K., Pont, F., Aigrain, S., et al. 2011, Monthly Notices of the Royal Astronomical Society, 416, 1443,  \ndoi: 10.1111/j.1365-2966.2011.19142.x  \nSkilling, J. 2004, 735, 395, doi: 10.1063/1.1835238  \n- -. 2006, Bayesian Analysis, 1, 833, doi: 10.1214/06-BA127 Soubiran, F., & Militzer, B. 2015, The Astrophysical Journal, 806, 228, doi: 10.1088/0004-637X/806/2/228\n- Teske, J. K., Wang, S., Wolfgang, A., et al. 2018, The Astronomical Journal, 155, 148, doi: 10.3847/1538-3881/aaab56\n- Thorngren, D. P., Fortney, J. J., Murray-Clay, R. A., & Lopez, E. D. 2016, The Astrophysical Journal, 831, 64, doi: 10.3847/0004-637X/831/1/64\n- Tovar Mendoza, G., Davenport, J. R. A., Agol, E., Jackman, J. A. G., & Hawley, S. L. 2022, AJ, 164, 17, doi: 10.3847/1538-3881/ac6fe6\n- Townsend, R., & Lopez, A. 2023, Journal of Open Source Software, 8, 4602, doi: 10.21105/joss.04602\n- Trotta, R. 2008, Contemporary Physics, 49, 71, doi: 10.1080/00107510802066753\n- Tsai, S.-M., Lee, E. K. H., Powell, D., et al. 2023, Nature, 617, 483, doi: 10.1038/s41586-023-05902-2\n- Venturini, J., Guilera, O. M., Haldemann, J., Ronco, M. P., & Mordasini, C. 2020, Astronomy & Astrophysics, 643, L1, doi: 10.1051/0004-6361/202039141  \nVenturini, J., Ronco, M. P., Guilera, O. M., et al. 2024, A fading radius valley towards M-dwarfs, a persistent density valley across stellar types, arXiv. http://arxiv.org/abs/2404.01967  \nWaalkes, W. C., Berta-Thompson, Z. K., Newton, E. R., et al. 2024, ApJ, 962, 97, doi: 10.3847/1538-4357/ad0bed Wallack, N. L., Batalha, N. E., Alderson, L., et al. 2024, JWST COMPASS: A NIRSpec/G395H Transmission Spectrum of the Sub-Neptune TOI-836c, arXiv. http://arxiv.org/abs/2404.01264  \nWelbanks, L., & Madhusudhan, N. 2019, AJ, 157, 206, doi: 10.3847/1538-3881/ab14de  \nWogan, N. F., Batalha, N. E., Zahnle, K., et al. 2024, JWST observations of K2-18b can be explained by a gas-rich mini-Neptune with no habitable surface, arXiv, doi: 10.48550/arXiv.2401.11082  \nYang, J., & Hu, R. 2024, Chemical mapping of temperate sub-Neptune atmospheres: Constraining the deep-interior H2O/H2 using the atmospheric CO2/CH4, arXiv. http://arxiv.org/abs/2406.01955  \n- Youngblood, A., France, K., Loyd, R. O. P., et al. 2016, ApJ, 824, 101, doi: 10.3847/0004-637X/824/2/101\n- Yurchenko, S. N., Mellor, T. M., Freedman, R. S., & Tennyson, J. 2020, Mon. Not. R. Astron. Soc., 496, 5282, doi: 10.1093/mnras/staa1874  \nZeng, L., Jacobsen, S. B., Sasselov, D. D., et al. 2019, Proceedings of the National Academy of Science, 116, 9723, doi: 10.1073/pnas.1812905116"}