Detecting atmospheres around planets with a radius below 1.6 R SUBSUB commonly referred to as rocky planets has proven to be challenging. However rocky planets orbiting M dwarfs are ideal candidates due to their favorable planettostar radius ratio. Here we present one transit observation of the SuperEarth L9859 d 1.58 R SUBSUB and 2.31 M SUBSUB at the limit of rockygasrich using the JWST NIRSpec G395H mode covering the 2.85.1 m wavelength range. The extracted transit spectrum from a single transit observation deviates from a flat line by 2.65.6 depending on the data reduction and retrieval setup. The hints of an atmospheric detection are driven by a large absorption feature between 3.3 and 4.8 m. A stellar contamination retrieval analysis rejected the source of this feature as being due to stellar inhomogeneities making the best fit an atmospheric model including sulfurbearing species suggesting that the atmosphere of L9859 d may not be at equilibrium. This result will need to be confirmed by the analysis of the second NIRSpec G395H visit in addition to the NIRISS SOSS transit observation.
['Exoplanet atmospheres', 'Exoplanet atmospheric composition', 'Transmission spectroscopy', 'Astronomy data reduction', 'Planetary atmospheres', 'Stellar atmospheres', 'Infrared spectroscopy', 'Super Earths', 'Extrasolar rocky planets', '487', '2021', '2133', '1861', '1244', '1584', '2285', '1655', '511', 'Astrophysics - Earth and Planetary Astrophysics']
Hints of a Sulfurrich Atmosphere around the 1.6 R SUBSUB SuperEarth L9859 d from JWST NIRspec G395H Transmission Spectroscopy
{'Hints of a sulfur-rich atmosphere around the 1.6 R ⊕ Super-Earth L98-59 d from JWST NIRSpec G395H transmission spectroscopy': "Am'elie Gressier , 1 N'estor Espinoza , 1, 2 Natalie H. Allen , 2 David K. Sing , 3, 4 Agnibha Banerjee , 5 Joanna K. Barstow , 5 Jeff A. Valenti , 1 Nikole K. Lewis , 6 Stephan M. Birkmann , 7 Ryan C. Challener , 8 Elena Manjavacas , 9, 3 Catarina Alves de Oliveira, , 7 Nicolas Crouzet , 10 and 1  \nTracy. L Beck  \n1 Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA  \n- 2 William H. Miller III Department of Physics and Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA\n- 3 Department of Physics and Astronomy, Johns Hopkins University, 3400 N. Charles Street, Baltimore, MD 21218, USA  \nDepartment of Earth and Planetary Sciences, Johns Hopkins University, 3400 N. Charles Street, Baltimore, MD 21218, USA  \n5 School of Physical Sciences, The Open University, Milton Keynes, MK7 6AA, UK  \n- 6 Department of Astronomy and Carl Sagan Institute, Cornell University, 122 Sciences Drive, Ithaca, NY 14853, USA  \n7 European Space Agency, European Space Astronomy Centre, Camino Bajo del Castillo s/n, E-28692 Villanueva de la Ca˜nada, Madrid, Spain  \n8 Department of Astronomy, Cornell University, 122 Sciences Drive, Ithaca, NY 14853, USA  \n9 AURA for the European Space Agency (ESA), ESA Office, Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA  \n10 Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands  \n(Accepted for publication in the Astrophysical Journal Letters, August 25, 2024.)", 'ABSTRACT': 'Detecting atmospheres around planets with a radius below 1.6 R ⊕ , commonly referred to as rocky planets (Rogers 2015; Rogers et al. 2021), has proven to be challenging. However, rocky planets orbiting M-dwarfs are ideal candidates due to their favorable planet-to-star radius ratio. Here, we present one transit observation of the Super-Earth L98-59 d (1.58 R ⊕ , 2.31 M ⊕ ), at the limit of rocky/gas-rich, using the JWST NIRSpec G395H mode covering the 2.8 to 5.1 µ m wavelength range. The extracted transit spectrum from a single transit observation deviates from a flat line by 2.6 to 5.6 σ , depending on the data reduction and retrieval setup. The hints of an atmospheric detection are driven by a large absorption feature between 3.3 to 4.8 µ m. A stellar contamination retrieval analysis rejected the source of this feature as being due to stellar inhomogeneities, making the best fit an atmospheric model including sulfur-bearing species, suggesting that the atmosphere of L98-59 d may not be at equilibrium. This result will need to be confirmed by the analysis of the second NIRSpec G395H visit in addition to the NIRISS SOSS transit observation.', '1. INTRODUCTION': "Pushing the instrumental limit to possibly detect an atmosphere around an Earth-like planet is the first step in understanding whether any of these planets could support life and to what extent we could detect it. The large planet-to-star radius ratio, coupled with the high frequency of the transit event, creates favorable conditions for searching for an atmosphere using transmission spectroscopy around planets orbiting M-dwarfs. However, in the search for an atmospheric signal around a rocky planet, one must balance the M-dwarf opportunity with the possibility that long exposure to X-ray and EUV radiations can strip away the atmosphere (Owen & Wu 2013; Owen 2019). Besides, M-dwarf stars can mimic  \nplanetary atmospheric features in the transmission spectrum from photospheric inhomogeneities, known as stellar contamination (Rackham et al. 2018).  \nSpace- and ground-based observations of rocky planets (R P < 1.6R ⊕ ) have not yet provided strong and compelling evidence of atmospheric signals. The Hubble Space Telescope Wide Field Camera 3 Grism 141 (HST WFC3 G141), widely used in transmission spectroscopy for its water band at 1.4 µ m, did not yield conclusive and/or significant evidence due to its precision and limited spectral range (de Wit et al. 2016, 2018; Wakeford et al. 2019; Mugnai et al. 2021; Edwards et al. 2021; Gressier et al. 2021; Garcia et al. 2022). Similarly, no strong constraints were found using ground-  \n4  \nbased (Diamond-Lowe et al. 2018) and Spitzer observations (Demory et al. 2011, 2012, 2015; Demory et al. 2016; Kreidberg et al. 2019; Zieba et al. 2022; Crossfield et al. 2022).  \nThe atmospheric features in transmission for warm rocky planets are expected to reach 20 to 100 ppm, which approaches the expected noise floor of the James Webb Space Telescope instruments(Greene et al. 2016). Even for the JWST, detecting an atmosphere around a terrestrial world is a challenge, and this has been confirmed within the first years of observations. To date, no strong evidence for an atmosphere around a rocky exoplanet has been found (Lustig-Yaeger et al. 2023; Moran et al. 2018; May et al. 2023; Kirk et al. 2024). TRAPPIST-1b and c photometric data points and transmission spectra are so far consistent with no or small atmospheric signals, with the interpretation of the transmission spectrum complicated by the presence of significant stellar contamination (Greene et al. 2023; Zieba et al. 2023; Lim et al. 2023).  \nThe L98-59 system is composed of 3 rocky transiting planets (Cloutier et al. 2019; Kostov et al. 2019), with radii below 1.58 R ⊕ , and a non-transiting planet (Demangeon et al. 2021) orbiting a bright M3-dwarf (K=7.01 mag). Spanning the HST and JWST cycles 1 and 2, the L98-59 system stands out as one of the most extensively studied multi-planetary systems in emission and transmission spectroscopy (cycle 1 GO and GTO programs 1201, 1224, 2512 and cycle 2 GO 3942, 3730, 4098). So far, the HST WFC3 G141 observations in the 1.1 to 1.7 µ m range have not provided strong evidence of an atmosphere for any planets in the system. Both Damiano et al. (2022) for planet b and Zhou et al. (2023) for planets c and d rejected cloud-free hydrogen and helium atmospheres. However, they could not rule out primary cloudy/hazy or water-rich atmospheres.  \nHere we present the result of the JWST-GTO-1224 program observations of the third planet of the system, L98-59 d, consisting of one transit covering the 2.87 to 5.17 µ m wavelength range. With a radius of 1.58 ± 0.08R ⊕ and a mass of 2.31 +0 . 46 -0 . 45 M ⊕ (Luque & Pall'e 2022), L98-59 d lies on the boundary defined by Rogers (2015) between rocky and gas-rich planets, making it a highly intriguing target. In this letter, we first present the observations in Section 2. In Section 3, we detail the two independent data analyses used to extract the transmission spectrum. Section 4 describes the interpretations we carried out, including atmospheric and stellar contamination forward and retrieval modeling. We discuss the results in the final Section 5.", '2. JWST OBSERVATIONS OF L98-59d': 'We observed one transit of L98-59 d using the NIRSpec/G395H mode as part of the JWST GTO cycle 1 Program 1224, led by PI Stephan Birkmann. The G395H instrument mode covers the wavelength range from 2.87 to 5.17 µ m, at a native resolution of R ∼ 2700. The grating is split over two detectors NRS1 and NRS2, with a gap between 3.72 to 3.82 µ m. The observations started on June 25 nd 2023, covered the full transit and sufficient baseline for a total of 5.34 hours, resulting in 2121 integrations, each composed of 6 groups up-theramp. The observation used the NIRSpec Bright Object Time Series mode with the NRSRAPID readout pattern, S1600A1 slit, and the SUB2048 subarray.', '3. DATA REDUCTION': 'The reduction of L98-59 d NIRSpec G395H observations was obtained using two independent pipelines, transitspectroscopy (Espinoza 2022) and FIREFLy (Rustamkulov et al. 2022, 2023) which have already been benchmarked for the JWST Early Release Science (ERS) program (Alderson et al. 2023; Rustamkulov et al. 2023). FIREFLy have provided data reductions for smaller planets observations with the JWST NIRSpec G395H (Lustig-Yaeger et al. 2023; Moran et al. 2023; May et al. 2023) which makes it a good comparison.', '3.1. Data reduction with transitspectroscopy': "We used the jwst pipeline version 1.12.5 and the transitspectroscopy pipeline (Espinoza 2022) 1 to reduce the NIRSpec G395H transit of L98-59d. The reduction process followed the jwst pipeline from uncal.fits to rateints.fits, excluding the jump step detection, which we customized using the transitspectroscopy algorithm. This algorithm identifies outliers at the pixel level within each group by calculating group differences, applying a median filter, and detecting jumps. Spectral tracing was performed using the trace spectrum function, smoothed with a spline function. We determined the trace position using a cross-correlation method with a Gaussian input, covering specified pixel ranges for NRS1 and NRS2. Background estimation involved masking pixels around the trace and calculating the median of the remaining pixels. To address 1/f noise (Jakobsen et al. 2022; Birkmann et al. 2022) - noise arising from the electronics of the readout pattern, which appears as column striping in the subarray image - we subtracted the median outof-transit frame and iteratively removed the median of  \nFigure 1. Top: Raw white light curves for NRS1 (left) and NRS2 (right) from transitspectroscopy and spectral light curves map. Middle : Corrected white light curves (colored data points) overplotted with the best-fit transit model (black line) from juliet . Bottom : residuals between raw light curves and transit model.  \n<!-- image -->  \nnon-omitted pixels near the trace position. The stellar spectrum was extracted using the getSimpleSpectrum routine with a 3-pixel radius aperture, minimizing outof-transit flux. We employed a box-extraction method to sum flux within the aperture and replaced outliers exceeding a 5 σ threshold with a 1D median filter. White light curves and pixel-level light curves were generated from the time series of 1D stellar spectra for both detectors.  \nWe independently fitted white and spectral light curves from each detector using the juliet Python package (Espinoza et al. 2019), employing nested sampling via dynesty (Speagle 2020). We restricted the time series observations to +1.0 hours from the start, as the transit begins in the first half-hour (see top panel Figure 1). The period ( P = 7 . 4507245 days), eccentricity ( e = 0), and argument of periastron ( ω = 90), are fixed based on values from (Luque & Pall'e 2022; Demangeon et al. 2021). We fitted the planet-to-star radius ratio, mid-transit time, impact parameter, and semi-major axis ratio, including limb-darkening coefficients determined with a square-root law and a uniform prior. Additionally, we included a mean-out-oftransit offset ( mflux ) and a jitter parameter( sigma w )  \nfor white noise. To account for correlated noise we used Gaussian Processes (GP) using the george package with a Mat'ern 3/2 kernel (Ambikasaran et al. 2015; ForemanMackey et al. 2017). This method captures complex systematics and has been validated in ERS WASP-39b studies (Feinstein et al. 2023; Alderson et al. 2023). A log-uniform prior was set for the GP amplitude between 0.001 and 1000 ppm, and an exponential prior was set for the GP length-scale, using time and the full-width half maximum (FWHM) as GP regressors. The white light curve fitting results are shown in Figure 1, and the detailed best-fit parameters and confidence intervals are provided in Appendix Table 2. We determined the optimal systematics detrending model based on Bayesian log-evidence, finding that the GP model with both time and FWHM as regressors offers the best fit. We evaluated several detrending approaches, including linear models with time and FWHM, and a GP with time only using the Mat'ern 3/2 kernel via the celerite package (Foreman-Mackey et al. 2017). For NRS1, the logevidence values were 6700 for the transit-only model, 6778 for the linear time model, 6780 for the linear time plus FWHM model, 6782 for the GP with time only, and 6784 for the GP with both time and FWHM. Similarly, for NRS2, the GP model including both time and  \nFigure 2. Transmission spectra of L98-59 d obtained using transitspectroscopy (pink) and FIREFLy (blue) pipelines. The spectral light curves from the two indepedent reductions were fitted using juliet with a similar parametrization at pixel resolution and then binned to R ∼ 100. The residuals between the two reductions show a global offset between the two reductions for the NRS2 detector.  \n<!-- image -->  \nFWHM achieved the highest log-evidence value of 6405, compared to 6302 for the transit-only model, 6212 for the linear time model, 6219 for the linear time plus FWHM model, and 6282 for the GP with time only.  \nFor the spectral light curves, we used a similar setup but fixed the mid-transit time, impact parameter, and semi-major axis ratio based on the best-fit results from the white light curves. The fitted values are T 0 (BJD TDB ) = 2460121 . 112518755 ± 7 × 10 -5 , b = 0 . 925 ± 0 . 004, and a/R ⋆ = 37 . 11 ± 0 . 55. Limbdarkening coefficients were determined using a truncated normal distribution centered on values estimated by the exotic (Grant & Wakeford 2022) package and PHOENIX models (Husser et al. 2013). Given the high impact parameter, constraining these coefficients was challenging, as discussed in Appendix B. The transmission spectrum was obtained by binning pixel-level transit depths to a resolution of R ∼ 100, resulting in 56 data points. Uncertainties in the final spectrum were calculated by determining the weighted mean and variance for  \neach bin, with errors derived from the weighted contributions of the high-resolution spectra's errors and their associated variance. This spectrum can be found in Appendix 3.", '3.2. Data reduction with FIREFLy': 'We reduced the L98-59 d data using the Fast InfraRed Exoplanet Fitting Lyghtcurve ( FIREFLy ; Rustamkulov et al. 2022, 2023) reduction suite. FIREFLy starts with a customised reduction using the STScI pipeline and the uncalibrated images and includes 1/ f destriping at the group level before the ramp is fit. We applied a custom superbias step, where the STScI pipeline superbias was scaled separately for NRS1 and NRS2 to match the trace-masked count level median value of each integration. This procedure was introduced in Moran et al. (2023) to help reduce offsets between the NRS1 and NRS2 detectors. The jump-step and dark-current stages of the STScI pipeline are skipped. We then use the custom-run 2D images after the gain scale step, and  \nperform customised cleaning of bad pixels, cosmic rays and hot pixels.  \nWe fit the transit light curves using a quadratic limbdarkening model ( q 1 , q 2 ; Kipping 2013, along with a polynomial function of time to remove overall baseline trends. In addition, we tested detrending the lightcurves with the X and Y positions on the detector and the scaling factor used in the superbias step. We found the NRS1 detector to have a step trend which evolved over the first 0.04 days requiring a polynomial in time with orders 1, 2 , 5 and 6. The higher-order trends modelled the trends seen at the beginning of the observation well. Conversely, NRS2 only required a second-order polynomial though a trend in Y-position was found to be needed. In addition a trend for both NRS1 and NRS2 with the superbias scaling factor was found to be needed. We first fit for both q 1 and q 2 , but found q 2 to be weakly constrained and consistent with 0, so q 2 was fixed to zero throughout the rest of the analysis. a/R ⋆ , T 0 , and the impact parameter b were fit for both NRS1 and NRS2, then fixed to the weighted-average value of both detectors. The spectroscopic fits followed the same procedure. However, q 1 was found not to have a strong wavelength dependence when fit spectroscopically and was weakly constrained, so it was fixed to the weighted average value for all wavelength channels ( q 1 =0.117 ± 0.017). The lightcurves were binned to produce a final transmission spectra with a resolution near R =60 with a total of 38 points.', '3.3. Data reductions comparison': "Wefit the pixel-light curves from transitspectroscopy and FIREFLy reductions with the parametrization described in Section 3.1. The results are in Figure 2. While both spectra exhibit similar shapes, a systematic offset is observed for NRS2 detector transit depths between the two reductions, with FIREFLy 's transmission spectrum exhibiting higher transit depth values than transitspectroscopy . Specifically, there is a higher median of residuals for NRS2 (+46ppm) compared to NRS1 ( -8ppm), potentially attributable to a different background treatment for NRS2 detector. The root mean square of the residuals is found to be 56 ppm for NRS2 and 20 ppm for NRS1. Similarly, the light curve fitting setup described in Section 3.2 is applied to transitspectroscopy 's light curves (see AppendixB.6).", '4. INTERPRETATION': 'We interpret the results using forward and retrieval modelling focusing on the transmission spectrum obtained with transitspectroscopy at R ∼ 100.', '4.1. Evidence of an atmosphere from 1D forward modelling': "We used Exo-REM (Baudino et al. 2015; Charnay et al. 2018; Blain et al. 2021) to create a grid of transmission spectra for L98-59 d. Exo-REM is a 1D selfconsistent radiative-convective code for exoplanets and brown dwarfs' atmospheres. The atmosphere is modelled using 80 layers between 10 -3 and 10 7 Pa, with an Eddy Diffusion coefficient of 10 -8 cm 2 /s -1 . We simulated a 300 × solar atmosphere, including 13 absorbing species in a hydrogen-rich atmosphere, considering both equilibrium and out-of-equilibrium chemistry. The metallicity is the factor by which all the elemental abundances except H are multiplied compared to their solar abundances (Lodders 2010). For out-of-equilibrium chemistry, Exo-REM incorporates disequilibrium processes for CH 4 -CO, CO-CO 2 , and N 2 -NH 3 -HCN reactions based on Zahnle & Marley (2014). The quenching level is determined by comparing the reaction timescale to the mixing time H 2 /K zz , where H is the atmospheric scale height and K zz the Eddy diffusion coefficient. The abundance of the species is governed by the thermochemical equilibrium below the quenching level, where the temperature and pressure are high enough so that the kinetics dominates. The internal temperature was set to 70 K, controlling the position of the convective layer. We also tested a pure H 2 O, CH 4 , CO 2 and H 2 S atmosphere.  \nFigure 3 shows the transmission spectrum obtained with transitspectroscopy and grid models from Exo-REM , along with corresponding χ 2 statistics. The pure H 2 S model provides the best fit to the spectrum with a chi-square difference of around 19 with the flatline. This represents a statistically significant improvement in fit, exceeding 3 σ (dof = 56). However, the observed variations in the transmission spectrum for a pure H 2 S model, which show differences of approximately 40 ppm, exceed the simple scale height prediction of ∼ 5 ppm. This is estimated using a 416 K equilibrium temperature, 34 g/cm 3 mean molecular weight. This discrepancy may arise from several factors including the modeling of the atmosphere with ExoREM, which simulates a pure H 2 S atmosphere within a hydrogendominated framework, potentially resulting in an inflated effective scale height. Additionally, the use of cross-section opacities made for hydrogen dominated atmosphere and detailed atmospheric physics in our models may explain the larger signal compared to the simplified scale height calculation. Further investigation into the support for sulfur-bearing species opacities in the atmospheric model will be conducted using a Bayesian framework in Section 4.2. Cloud-free models with 300  \nFigure 3. Transmission spectrum of L98-59 d obtained using transitspectroscopy and binned to R ∼ 100 (black points) compared to forward models from Exo-REM (coloured dashed and filled lines). The pure H 2 S model provides the best-fit according to χ 2 statistics. No offset is added between NRS1 and NRS2.  \n<!-- image -->  \ntimes solar abundance exhibit large CH 4 features at 3.3 µ m and CO/CO 2 features around 4.5 µ m, which are not observed in the L98-59 d transmission spectrum. Despite the pure H 2 S model being a likely non-physical atmospheric model, it is the best-fit according to χ 2 statistics. The modeled flat-line presented in Figure 3 is computed using the results of the retrieval analysis with TauREx for consistency (see Section 4.2). Using a weighted mean of spectrum points yields a chi-squared of 80.47, instead of 81.62. This results might indicate an atmosphere with the presence of sulfur-bearing species or stellar contamination.", '4.2.1. Retrieval configurations': "We conduct an atmospheric and a stellar contamination retrievals of L98-59 d's properties based on its observed spectra obtained through transitspectroscopy using TauREx 3 (Al-Refaie et al. 2021) 2 and exoretrievals (Espinoza et al. 2019).  \nFor the atmospheric retrieval analysis, we used the Multinest algorithm (Feroz et al. 2009; Buchner et al.  \n2014) with 100 layers in the atmosphere, spanning from 10 -3 to 10 6 Pa. The model included an evidence tolerance of 0.5 and 1500 live points. Stellar parameters were fixed (radius: 0.303 R ⊙ , T star : 3415 K, metallicity: -0.46). The planetary radius was fitted within the range of 0.79 to 2.37 R ⊕ to explore the full parameter space, account for uncertainties including limb-darkening effects, and ensure robustness in the model. The temperature-pressure profile ranged from 200 to 600 K. Clouds were modeled as grey clouds with a top pressure between 10 -3 and 10 6 Pa. The helium-to-hydrogen ratio was set to 0.17, and an offset between NRS1 and NRS2 transmission spectra was fitted uniformly between -100 and +100 ppm using the taurex-offset plugin 3 .  \nEquilibrium chemistry : We used ggchem 4 (Woitke et al. 2018) for equilibrium chemistry, including molecular line lists and continuum from the ExoMol project (Tennyson et al. 2016; Chubb et al. 2021), HITEMP (Tennyson & Yurchenko 2018), and HITRAN (Rothman et al. 1987; Rothman et al. 2010). The active molecules included are : H 2 O (Polyansky et al. 2018), CO (Li et al. 2015), CO 2 (Rothman et al.  \n2010), CH 4 (Yurchenko et al. 2017), NH 3 (Yurchenko et al. 2011), H 2 S (Azzam et al. 2016), HCN(Barber et al. 2014), C 2 H 2 (Chubb et al. 2020), C 2 H 4 (Mant et al. 2018), SiO(Yurchenko et al. 2021), SO 2 (Underwood et al. 2016), TiO(McKemmish et al. 2019), and VO(McKemmish et al. 2016). Collision-induced absorption (CIA) and Rayleigh scattering were also considered. We fitted for C/O and S/O ratios (0.1 to 2) and atmospheric metallicity (1 to 1000 times solar).  \nFree chemistry : The active molecules included are H 2 O, CO, CO 2 , CH 4 , SiO, H 2 S, HCN, NH 3 , and SO 2 . CIA and Rayleigh opacities were included. We fitted the abundance of each molecule (log(X VMR )) between 10 -12 and 0.3 (between -12 and -0.5 in log space), ensuring the total molecular abundances never exceed one.  \nFree chemistry with N 2 : The setup is similar to the free chemistry retrieval but included N 2 as an inactive gas. The N 2 /H 2 ratio was fitted between 0.001 and 2.  \nFlat-line : The flat-line is a retrieval with no opacity sources, fitting for planetary radius and temperature. This was used to assess the significance of the atmospheric detection using Bayesian evidence, evaluated using the Bayes's theorem. The Bayes factor, positively designed is computed between an atmospheric model and the flat-line and then converted in σ significance using the formalism of Trotta (2008) and Benneke & Seager (2013). This retrieval is performed for each atmospheric retrieval on the transmission spectrum with the retrieved offset applied and fixed.  \n4.2.2. Retrieval results : an atmosphere around L98-59d ?  \nThe retrieval results using transitspectroscopy reduction are summarised in Table 1. We detected an atmosphere around L98-59 d at 5.6 σ with the free chemistry retrieval and at 2.7 σ with the equilibrium chemistry retrieval. Both retrievals report a high abundance of H 2 S: log(H 2 S)= -0 . 74 +0 . 14 -0 . 49 (free chemistry) and log(H 2 S)= -0 . 21 +0 . 07 -0 . 09 (equilibrium), with Bayesian evidences log Z = 464 . 99 and log Z = 460 . 31, respectively. The flat-line fits yield log Z = 451 . 19 and log Z = 458 . 00.  \nThe free chemistry retrieval best-fit model is presented in Figure 4 and the corresponding posterior distributions are in Appendix C.7 while the equilibrium chemistry retrieval is presented in Appendix C.8. The results align with Paper II using NEMESISPY(Irwin et al. 2008; Yang et al. 2023). An offset of 37 . 6 +15 . 0 -16 . 4 ppm  \nwas found between NRS2 and NRS1 with transit depths higher for NRS2. Both models suggest significant H 2 S opacity, with CO 2 (equilibrium) and SO 2 (free chemistry) as additional opacity contributers. The free chemistry model fits significantly better than the equilibrium model by 3.8 σ suggesting that L98-59 d's atmosphere is not at equilibrium. The SO 2 VMR in the free retrieval shows a bi-modal distribution, with a constraint either above -1 . 0 or an upper limit below -1 . 5. However, we can put upper limits on every other molecules of the fit. In particular, H 2 O's VMR is constrained below 10 -5 suggesting a water-poor atmosphere for this planet that has density compatible with a water-rich composition (Luque & Pall'e 2022). Excluding both H 2 S and SO 2 from the free chemistry retrieval indicates that the model with these molecules is favored at 4.5 σ using model comparison with the Bayes factor.  \nAdding N 2 as an inactive gas in the free chemistry retrieval did not significantly alter results, favoring the atmospheric model over a flat-line at 4.5 σ . The only strongly constrained parameters are the planetary radius at 10 bars, the log VMR of H 2 S, toward the edge of the prior (0.3), and the top pressure of the grey clouds. The log VMR of SO 2 is tentatively constrained. For this reason, we run a fit with H 2 S and SO 2 as the only active molecular opacity sources, increasing the prior range to 1, and fit for their molecular abundances while fixing the temperature to the equilibrium temperature at 416K. This atmospheric model is preferred at 5.3 σ to a flat-line.  \n4.2.3. Stellar contamination : is there an atmosphere or stellar inhomogeneities ?  \nOne of the early results from observations of small planets with JWST is how important it is to consider heterogeneities on the surface of the host star when looking for atmospheric signals in transmission. Multiple small planet observations have shown contamination from the presence of hot and cold spots on the photosphere of the host star during transit, which imprint false slopes and can potentially even inject false molecular features into the transmission spectrum (see e.g. Moran et al. 2023; Lim et al. 2023). This is especially important for planets orbiting small M stars like L9859, which are known to show a high level of magnetic activity and associated surface features (Rackham et al. 2018).  \nTo test if the spectrum could be explained by stellar contamination caused by inhomogeneities (hot and/or cold spots) on the stellar surface, we carry out retrievals using exoretrievals (Espinoza et al. 2019). Based on the framework presented in Rackham et al. (2018), this allows us to model the spectrum created by occulted  \nFigure 4. Top panel: Transmission spectrum of L98-59 d obtained using transitspectroscopy and binned to R ∼ 100 (black points) compared to a free chemistry retrieval model from TauREx (purple) and NEMESISPY (dashed blue) from Banerjee et al. 2024, hereafter Paper II. Middle : Opacity contributions from the best-fit model. The model suggests a contribution from H 2 S (purple) and SO 2 (green) opacities. Bottom : Transmission spectrum compared to stellar (yellow) and atmospheric (purple) retrievals models.  \n<!-- image -->  \nTable 1. Summary of Atmospheric retrieval results using TauREx on L98-59d transmission spectrum obtained with transitspectroscopy 's reduction at R ∼ 100.  \nand/or unocculted spots on the stellar surface during the transit, by replacing a portion of the surface with a stellar model of another temperature. Since we do not see any evidence of spot crossing events in the light curves, we only consider unocculted spots in our retrieval. We consider the cases of only hot spots, only cold spots, or both being present on the stellar surface, and allow these spot temperatures to vary over a large parameter space, from 1800 K to 3250 K and from 3550 K to 5000 K, respectively, for a photospheric temperature of 3415 ± 135. We also place a wide prior on the spot covering fraction, allowing the spots to together cover up to the entirety of the unocculted photosphere. We also tested using both default PHOENIX models and BT-SETTL models considering the relatively low photospheric temperature of L98-59, but saw no appreciable difference in the results. Even with this very large explored parameter space, none of the stellar contamination models are able to explain the features seen in the spectrum, especially those at the longest wavelengths. Figure 4 bottom panel shows the best-fit model from the stellar contamination retrieval results. The corresponding posterior distributions are in Appendix C.10. We also show the best-fit model from an atmospheric retrieval analysis with TauREx . In these fits, no offset between NRS1 and NRS2 is applied. The H 2 S/SO 2 rich atmosphere is favored compared to the stellar contamination suggesting the large feature between 3.3 to 4.8 µ m, if confirmed by follow-up observations is originated from the atmosphere.", '4.2.4. Comparison between data reductions': "To confirm the atmospheric detection, we applied the free chemistry retrieval setup to the FIREFLy transmission spectrum at R ∼ 100 and to the spectrum obtained at R ∼ 60 from an independent fitting method (see Figure 2). The results of these fits are shown in Appendix Figure C.11. The atmospheric detection was less significant with FIREFLy , yielding 3.0 σ for R ∼ 100 and 2.6 σ for R ∼ 60 spectra.  \nThe discrepancy at R ∼ 100 primarily arises from differences at longer wavelengths in the NRS2 detectors. transitspectroscopy data showed a decrease in transit depth between 4.5 and 5 µ m, fitted with H 2 S opacity (see the middle panel in Figure 4) while FIREFLy 's transmission spectrum remains flat across the entire wavelength range of the NRS2 detector. However, sulfurbearing species consistently appeared as the primary source of opacity in all retrievals. The NRS1 and NRS2 offset varied by reduction method: +37 ppm for transitspectroscopy and -14 ppm for FIREFLy , both consistent with zero at 1 σ . Removing this offset is un-  \nly to significantly change the FIREFLy result given the small offset within 1 σ uncertainty.  \nThe R ∼ 60 spectrum from FIREFLy with an independent light curve fitting showed a tentative atmospheric detection (2.6 σ ). This finding raises concerns, suggesting that the detection made with other binning and light curve fitting may be attributed to random noise. Specifically, the discrepancy observed at shorter wavelengths, between 3 and 3.3 µ m, between the two reductions could account for this result (see Figure B.6). The variability observed in FIREFLy 's transmission spectrum in this wavelength range resulted in a seemingly flatter spectrum and a statistically non-detection, whereas transitspectroscopy 's data indicated a pronounced decrease in transit depths, well-fitted by H 2 S's opacity (see the middle panel in Figure 4). We note that employing the same light curve fitting and binning, the decrease is observed in FIREFLy 's reduction as well (see Figure 2).  \nIn AppendixC, we assess the impact of introducing an offset between NRS1 and NRS2 transit depths on atmospheric detection for L98-59 d (see Appendix C.5) and evaluate the influence of spectral resolution (see Appendix C.6).", '5. CONCLUSIONS': "We presented one transit observation of the SuperEarth L98-59 d with the JWST NIRSpec G395H mode. Our study highlights several challenges encountered in the data reduction process when analyzing transmission spectra of small planets. First, the limb-darkening coefficients pose a challenge for grazing planets due to their impact on the light curve fitting process. Even though our two data reductions agree when using the same limb-darkening coefficients, remaining discrepancies at specific wavelength led to different interpretation in the retrieval analysis.  \nFrom this single transit, our analysis suggests no stellar contamination but hints at potential atmospheric detection, though significance varies with data reduction methods. In the best case scenario, we detected an atmosphere with a 5.6 σ significance. An independent analysis led to a 2.6 σ tentative detection. Atmospheric retrieval models suggest the presence of sulfur-bearing species with hydrogen and helium as background gases, although the inferred high abundance levels are not yet well understood. The discovery of sulfur dioxide in the atmospheres of hot Jupiters (Rustamkulov et al. 2023; Alderson et al. 2023; Powell et al. 2024) and Neptunelike planets(Dyrek et al. 2023; Benneke et al. 2024; Holmberg & Madhusudhan 2024) are key-result in early JWST analysis being evidence of photo-chemistry. Sul-  \nr species and clouds are expected in giant planets and brown-dwarf atmosphere (Morley et al. 2012; Tsai et al. 2023; Polman et al. 2023). Crossfield (2023) showed that the ability to detect SO 2 in exoplanet atmospheres provides a crucial test for different planet formation models, revealing that sulfur's volatility and abundance can distinguish between planetesimal and pebble-accretion scenarios. Photochemistry could also produce H 2 S and SO 2 in terrestrial atmosphere as it was predicted by Hu et al. (2013). The presence of sulfur-bearing species in rocky planets could also be explained by out-gassing, volcanism, interaction between the atmosphere and the rocky surface. The recent study of Janssen et al. (2023) investigated sulfur's presence and detectability in rocky exoplanet atmospheres using thermochemical equilibrium models at the crust-atmosphere interface, considering surface temperatures of 500-5000 K and pressures of 1-100 bar, with various element abundances based on common rock compositions. They showed that at temperatures between 1000 and 2000 K, gaseous sulfur concentrations can reach up to 25%. They conclude that the most abundant sulfur molecules are SO, SO 2 , H 2 S, and S 2 with potentially detectable features in transmission spectra around 4 µ m, 7-8 µ m, and beyond 15 µ m.  \nThe difference between the data reductions and the uncertainties on the data points are magnified when dealing with observations from a single visit. (May et al. 2023) highlighted for GJ 1132 b the importance of multiple visits when claiming atmospheric detection for rocky planets. Our transmission spectrum might have random noise fluctuations that suggested the detection of an atmosphere which emphasizes the need for multiple transits observations. This result will have to be compared to the second G395H visit from Cycle 2 GO program 4098 (PI: Benneke) but also to the NIRISS SOSS transit from Cycle 1 GTO program 1201 (PI: Lafreni'ere).  \nOnly hints of atmospheric detections have been found so far for rocky planets with radii below 1.6 R ⊕ around M-dwarfs. Notably, planets from the GO program 1981 (PIs: Stevenson and Lustig-Yaeger), in order of increasing radius, include GJ 341 b with a radius of 0.92 R ⊕ (Kirk et al. 2024), LHS 475 b with 0.99 R ⊕ (LustigYaeger et al. 2023), GJ 1132 b with 1.1 R ⊕ (May et al. 2023), and GJ 486 b with 1.3 R ⊕ (Moran et al. 2023), yet they have not provided compelling evidence of an atmosphere. However, similar JWST observations on planets with radii above 1.6R ⊕ , such as K2-18 b (2.61 R ⊕ ) (Madhusudhan et al. 2023) and TOI-270d (2.31R ⊕ ), 55 Cancri e (1.97 R ⊕ )(Hu et al. 2024) have provided strong molecular detections. If confirmed, the detection of sulfur-bearing species in an hydrogen-dominated atmosphere around L98-59 d, a planet with a radius of 1.58 Earth radii, would be a significant result, as it lies right at the cutoff predicted by Rogers (2015); Rogers et al. (2021) for planets to have retained their primary hydrogen-helium atmosphere. Planets with a radius above 1.6R ⊕ have a density too low to be compatible with a silicate and iron-only core composition. A second study Paper II, from the NIRSpec GTO collaboration will explore in more details the possible atmospheric scenarios for L98-59 d.  \nN.H.A. acknowledges support by the National Science Foundation Graduate Research Fellowship under Grant No. DGE1746891.  \nThe JWST data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute. The specific observations analyzed can be accessed via 10.17909/nrxs-cx46.  \nFacilities: JWST NIRSPEC G395H", 'REFERENCES': "Al-Refaie, A. F., Changeat, Q., Waldmann, I. P., & Tinetti, G. 2021, The Astrophysical Journal, 917, 37,  \ndoi: 10.3847/1538-4357/ac0252 Alderson, L., Wakeford, H. R., Alam, M. K., et al. 2023, Nature, 614, 664-669, doi: 10.1038/s41586-022-05591-3 Ambikasaran, S., Foreman-Mackey, D., Greengard, L., Hogg, D. W., & O'Neil, M. 2015, IEEE Transactions on Pattern Analysis and Machine Intelligence, 38, 252, doi: 10.1109/TPAMI.2015.2448083  \nAzzam, A. A. A., Tennyson, J., Yurchenko, S. N., & Naumenko, O. V. 2016, Monthly Notices of the Royal Astronomical Society, 460, 4063-4074,  \ndoi: 10.1093/mnras/stw1133  \nBarber, R. J., Strange, J. K., Hill, C., et al. 2014, Monthly Notices of the Royal Astronomical Society, 437, 1828,  \ndoi: 10.1093/mnras/stt2011 Baudino, J.-L., B'ezard, B., Boccaletti, A., et al. 2015, Astronomy & Astrophysics, 582, A83, doi: 10.1051/0004-6361/201526332 Benneke, B., & Seager, S. 2013, The Astrophysical Journal, 778, 153, doi: 10.1088/0004-637x/778/2/153 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. https://arxiv.org/abs/2403.03325  \nUnderwood, D. S., Tennyson, J., Yurchenko, S. N., et al.  \n2016, Monthly Notices of the Royal Astronomical  \nSociety, 459, 3890-3899, doi: 10.1093/mnras/stw849  \nWakeford, H. R., Lewis, N. K., Fowler, J., et al. 2019, The Astronomical Journal, 157, 11,  \ndoi: 10.3847/1538-3881/aaf04d  \n- Woitke, P., Helling, C., Hunter, G. H., et al. 2018, Astronomy & Astrophysics, 614, A1,  \ndoi: 10.1051/0004-6361/201732193  \n- Yang, J., Irwin, P. G. J., & Barstow, J. K. 2023, Testing 2D temperature models in Bayesian retrievals of atmospheric properties from hot Jupiter phase curves. https://arxiv.org/abs/2305.10249  \nYurchenko, S. N., Amundsen, D. S., Tennyson, J., & Waldmann, I. P. 2017, Astronomy & Astrophysics, 605, A95, doi: 10.1051/0004-6361/201731026  \n- Yurchenko, S. N., Barber, R. J., & Tennyson, J. 2011, Monthly Notices of the Royal Astronomical Society, 413, 1828-1834, doi: 10.1111/j.1365-2966.2011.18261.x  \nYurchenko, S. N., Tennyson, J., Syme, A.-M., et al. 2021, Monthly Notices of the Royal Astronomical Society, 510, 903-919, doi: 10.1093/mnras/stab3267  \nZahnle, K. J., & Marley, M. S. 2014, The Astrophysical Journal, 797, 41, doi: 10.1088/0004-637x/797/1/41 Research in Astronomy and Astrophysics, 23, 025011,  \n- Zhou, L., Ma, B., Wang, Y.-H., & Zhu, Y.-N. 2023, doi: 10.1088/1674-4527/acaceb  \nZieba, S., Zilinskas, M., Kreidberg, L., et al. 2022, Astronomy & Astrophysics, 664, A79, doi: 10.1051/0004-6361/202142912  \nZieba, S., Kreidberg, L., Ducrot, E., et al. 2023, Nature, 620, 746-749, doi: 10.1038/s41586-023-06232-z", 'A. DATA REDUCTION : WHITE LIGHT CURVE BEST-FIT RESULTS': 'Table 2. L98-59 d White Light Curve fitting Results for the transitspectroscopy reduction using juliet .  \nDescription of parameters:Mid-transit time (BJD TDB ). Ratio of the semi-major axis to the stellar radius. Planet to star radius ratio.  \nFirst limb-darkening coefficient of the squareroot law with the Kipping (2013) parametrization. Second limb-darkening coefficient. Out-of-transit flux offset. White noise jitter parameter. GP amplitude. GP length scales.', 'B. DATA REDUCTION: TREATMENT OF THE LIMB-DARKENING COEFFICIENTS': "L98-59 d had a high impact parameter b > 0 . 9, which makes it challenging to constrain the limb-darkening coefficient in the light curve fitting. The choice of the limb-darkening coefficient treatment impact the extracted transit spectrum particularly the absolute transit depth value. First, we ran tests: one involving fitting the limb-darkening coefficients uniformly between 0 and 1 using a quadratic law, and another using a square-root law. For this test, we fit the light curves from the transitspectroscopy 's reduction at pixel resolution and bin the spectrum after. Results of this test are in Appendix Figure B.5. We note that there is a 400ppm offset between the two resulting transit spectra. The fit involving the quadratic law has a higher mean transit depth and lower precision. The shape of the spectrum is the same but the transit depth values and precision are sensitive to the choice of the limb-darkening treatment.  \nAdditionally, we perform a fit with fixed limb-darkening coefficients ( q 1 = 0 . 13 and q 2 = 0 . 0), derived from the combined best-fit of NRS1 and NRS2 white light curves, with q 2 = 0 . 0 being unconstrained. This parametrization is applied to pixel-level transitspectroscopy light curves, while FIREFLy light curves are binned to a R ∼ 60 and fitted with these fixed coefficients and an independent detrending linear model. This fitting is described in Section 3.2. The results of this additionnal test is in Appendix B.6.  \nFigure B.5. Top: Transmission spectra of L98-59 d obtained using the reductions from transitspectroscopy and binned to a resolution of R ∼ 60. The two spectra are derived by fitting the limb-darkening coefficients uniformly between 0 and 1 using a quadratic law (orange) and employing a square-root law (blue). The light orange points are obtained by applying an offset corresponding to the median differences between the two spectra. Bottom: Residuals between the two transmission spectra after subtraction of the median transit depth. Each detector is treated independently. The shape of the two spectra is similar but there is a 400ppm offset between the two reductions.  \n<!-- image -->", 'C. ADDITIONAL RETRIEVAL ANALYSIS AND FIGURES': "C.1. Free chemistry retrieval analysis  \n- C.2. Equilibrium chemistry retrieval analysis\n- C.3. Stellar contamination retrieval analysis  \nC.4. Atmospheric retrieval analysis on the FIREFLy data reduction  \nC.5. Does the offset between NRS1 and NRS2 impact the atmospheric detection ?  \nTo assess the influence of introducing an offset between the NRS1 and NRS2 transit depth values, we conducted an atmospheric retrieval fit on the transmission spectrum of L98-59 d at R ∼ 100 obtained from transitspectroscopy . The best-fit model is shown in Figure 4 bottom panel along with the stellar retrieval model. The atmospheric model is preferred at 3.4 σ compared to a flat-line with SO 2 and H 2 S as the main opacities. Moreover, incorporating these two molecules is favored at a significance of 3.2 σ compared to a model without them. Introducing the offset led to higher NRS2 transit depths in transitspectroscopy 's data, thereby reinforcing the atmospheric detection, particularly of SO 2 at 4.2 µ m. However, removing the offset did not impact the detection of an atmosphere in transitspectroscopy 's reduction.  \nFigure B.6. Transmission spectra of L98-59 d obtained using transitspectroscopy (pink) and FIREFLy (blue) pipelines. The spectral light curves from the two independent reductions were fitted independently. transitspectroscopy are fitted at pixel-level resolution using Gaussian Processes to model the noise and then binned to the native resolution of the FIREFLy 's reduction to approximately a R ∼ 60. FIREFLy light curves are binned to R ∼ 60 and then fitted using a linear de-trending model. For both fitting the coefficients of the limb-darkening coefficients are fixed to q 1 = 0 . 13 and q 2 = 0 . 0.  \n<!-- image -->", 'C.6. Does the spectral resolution impact the atmospheric detection ?': "To examine the impact of spectral resolution on our results, we applied the atmospheric free chemistry retrieval setup to transitspectroscopy 's data binned at different resolutions, following the fit at pixel-level resolution. We then evaluated the atmospheric detection by comparing Bayesian evidences with a flat-line, representing a fit without opacity sources. We did not fit for an offset between NRS1 and NRS2 in this analysis. The choice of spectral resolution did not significantly affect the atmospheric detection obtained with transitspectroscopy 's reduction.  \nFigure C.7. Posterior distributions from the free retrieval analysis. Not all parameters in the fit are plotted for clarity.  \n<!-- image -->  \nTable 3. L98-59 d extracted transit depths from transitspectroscopy at R ∼ 100  \nFigure C.8. Top panel: Transmission spectrum of L98-59 d obtained using transitspectroscopy and binned to R ∼ 100 (black points) compared to an equilibrium chemistry retrieval model from TauREx (light green). Bottom : Opacity contributions from the best-fit model. The model suggests a contribution from H 2 S (purple) and CO 2 (red) opacities.  \n<!-- image -->  \nFigure C.9. Posterior distributions from the equilibrium chemistry retrieval analysis. Not all parameters in the fit are plotted for clarity.  \n<!-- image -->  \nFigure C.10. Posterior distributions from the stellar retrieval analysis.  \n<!-- image -->  \nFigure C.11. Top Left:Transmission spectrum of L98-59 d obtained using transitspectroscopy (pink data points) and FIREFLy (blue data points) and binned to R ∼ 100 compared to atmospheric retrievals models with 1 σ uncertainty. An Offset is fitted between NRS1 and NRS2. Top Right : Posterior distributions from the atmospheric retrieval analyses. We do not represent all parameters for clarity. The values are printed for the FIREFLy atmospheric fit. Bottom Left :Transmission spectrum of L98-59 d obtained using FIREFLy (blue data points) binned to R ∼ 60 compared to an atmospheric retrieval models with 1 σ uncertainty. An Offset is fitted between NRS1 and NRS2. Bottom Right : Posterior distributions from the atmospheric retrieval analysis.  \n<!-- image -->"}