We report observations of the atmospheric transmission spectrum of the subNeptune exoplanet GJ 3470 b taken using the NearInfrared Camera on JWST. Combined with two archival Hubble Space TelescopeWideField Camera 3 transit observations and 15 archival Spitzer transit observations we detect water methane sulfur dioxide and carbon dioxide in the atmosphere of GJ 3470 b each with a significance of gt3. GJ 3470 b is the lowestmassand coldestexoplanet known to show a substantial sulfur dioxide feature in its spectrum at M SUB p SUB  11.2 M SUBSUB and T SUBeqSUB  600 K. This indicates that disequilibrium photochemistry drives sulfur dioxide production in exoplanet atmospheres over a wider range of masses and temperatures than has been reported or expected. The water carbon dioxide and sulfur dioxide abundances we measure indicate an atmospheric metallicity of approximately 100 solar. We see further evidence for disequilibrium chemistry in our inferred methane abundance which is significantly lower than expected from equilibrium models consistent with our measured water and carbon dioxide abundances.
['Exoplanet atmospheres', 'Exoplanet atmospheric composition', 'Exoplanet formation', '487', '2021', '492', 'Astrophysics - Earth and Planetary Astrophysics']
Sulfur Dioxide and Other Molecular Species in the Atmosphere of the SubNeptune GJ 3470 b
{'Sulfur Dioxide and Other Molecular Species in the Atmosphere of the Sub-Neptune GJ 3470 b': "Thomas G. Beatty, 1 Luis Welbanks, 2 Everett Schlawin, 3 Taylor J. Bell, 4 Michael R. Line, 2 Matthew Murphy, 3 Isaac Edelman, 4, 5 Thomas P. Greene, 6 Jonathan J. Fortney, 7 Gregory W. Henry, 8 Sagnick Mukherjee, 9 Kazumasa Ohno, 10, 9 Vivien Parmentier, 11 Emily Rauscher, 12 Lindsey S. Wiser, 2 and Kenneth E. Arnold 1  \n1 Department of Astronomy, University of Wisconsin-Madison, Madison, WI, USA  \n2 School of Earth and Space Exploration, Arizona State University, Tempe, AZ, USA  \n3 Steward Observatory, University of Arizona, Tucson, AZ, USA  \n4 Bay Area Environmental Research Institute, NASA's Ames Research Center, Moffett Field, CA, USA  \n5 Space Science and Astrobiology Division, NASA's Ames Research Center, Moffett Field, CA, USA  \n6 Space Science and Astrobiology Division, NASA's Ames Research Center, Moffett Field, CA, USA  \n7 Department of Astronomy and Astrophysics, University of California, Santa Cruz, 1156 High Street, Santa Cruz, CA 95064, USA 8 Tennessee State University, Nashville, TN, USA (retired)  \n9 Department of Astronomy and Astrophysics, University of California Santa Cruz, Santa Cruz, CA, USA  \n10 Division of Science, National Astronomical Observatory of Japan, Tokyo, Japan  \n11 Laboratoire Lagrange, Observatoire de la Cˆote d'Azur, Universit´e Cˆote d'Azur, Nice, France  \n12 Department of Astronomy, University of Michigan, Ann Arbor, MI, USA", 'ABSTRACT': 'We report observations of the atmospheric transmission spectrum of the sub-Neptune exoplanet GJ 3470 b taken using the Near-Infrared Camera (NIRCam) on JWST. Combined with two archival HST/WFC3 transit observations and fifteen archival Spitzer transit observations, we detect water, methane, sulfur dioxide, and carbon dioxide in the atmosphere of GJ 3470 b, each with a significance of > 3 σ . GJ 3470 b is the lowest mass - and coldest - exoplanet known to show a substantial sulfur dioxide feature in its spectrum, at M p =11.2 M ⊕ and T eq =600K. This indicates disequilibrium photochemistry drives sulfur dioxide production in exoplanet atmospheres over a wider range of masses and temperatures than has been reported (Tsai et al. 2023) or expected (Polman et al. 2023). The water, carbon dioxide, and sulfur dioxide abundances we measure indicate an atmospheric metallicity of approximately 100 × Solar. We see further evidence for disequilibrium chemistry in our inferred methane abundance, which is significantly lower than expected from equilibrium models consistent with our measured water and carbon dioxide abundances.', '1. INTRODUCTION': "One of the surprising results from exoplanet demographics studies is the prevalence of Neptune- and subNeptune-sized exoplanets (Fulton & Petigura 2018). Given their ubiquity, understanding how these exoplanets form is an important part of understanding planet formation. One way to do this is by studying their atmospheric compositions, since the carbon, oxygen, and overall metal abundances in Neptune- and subNeptune-mass planetary atmospheres are expected to be strongly influenced by their formation histories (Mordasini et al. 2016). Existing observations of these atmospheres from HST provide some water abundance measurements (Wakeford et al. 2017; Morley et al. 2017; Benneke et al. 2019a), but the difficulty in making these measurements and the lack of detections of carbonbearing molecules has made comparing to planet for-  \nmation theories difficult (Welbanks et al. 2019). Recent JWST observations of K2-18b (Madhusudhan et al. 2023) and TOI-270d (Benneke et al. 2024; Holmberg & Madhusudhan 2024) have begun to reveal spectra with strong molecular features in this regime, though many other sub-Neptunes have appeared featureless (LustigYaeger et al. 2023; May et al. 2023; Kirk et al. 2024).  \nOne of the canonical sub-Neptunes in this region of parameter space that is known to have features in its transmission spectrum is GJ 3470 b, which has a radius of 4.2 R ⊕ , a mass of 11.2 M ⊕ , and an equilibrium temperature of 600 K (Bonfils et al. 2012). This makes GJ 3470 b a relatively cool and relatively low-density member of the sub-Neptune population. Ground-based transit observations identified a Rayleigh-scattering slope in the atmosphere (Crossfield et al. 2013; Dragomir et al. 2015) and indicated that the atmosphere was likely methane-  \neted (Crossfield et al. 2013). HST and Spitzer observations of the transmission spectrum of GJ 3470 b (Benneke et al. 2019b) have shown significant water absorption near 1.4 µ m. These observations again indicated that the planetary atmosphere is likely depleted in methane, and the analysis also suggested GJ 3470 b has an atmospheric metallicity close to that of the Sun (as measured via the water abundance). These HST and Spitzer observations also suggested a possibly strong slope in the transmission spectrum between 2.5 µ m and 3.5 µ m, which was hypothesized to be caused by Mie scattering in the atmosphere.  \nThe apparent low metallicity of GJ 3470 b has been a puzzle, since both Uranus and Neptune have wellmeasured metallicities at approximately 100 times the Solar value (Atreya et al. 2016). Planet formation theories also suggest that the metallicity of a lower-mass planet like GJ 3470 b should be 200 to 300 times Solar (Mordasini et al. 2016). However, it is important to note that estimates for the metallicities of Solar System planets are based on methane-to-hydrogen ratios, while most previous metallicity estimates for exoplanets (including GJ 3470 b) have primarily used water-to-hydrogen ratios (Welbanks et al. 2019). This difference in metallicity estimation methods has made detailed comparisons between Solar System and exoplanet results difficult and has been necessitated by the lack of spectroscopic coverage of strong, isolated, methane absorption features in existing exoplanet atmosphere observations from either the ground or space.  \nGJ 3470 b is also an interesting target for atmospheric characterization observations because it appears to be significantly less cloudy than similar-mass exoplanets orbiting bright stars. Specifically, the other two 'canonical' exoplanets in this mass range, GJ 436 b and GJ 1214 b, both show infrared transmission spectra that are consistent with flat lines in both HST/WFC3 (Knutson et al. 2014; Kreidberg et al. 2014) and more recent JWST/MIRI (in the case of GJ 1214 b, Gao et al. 2023) observations, indicating that both have high-altitude aerosols around their terminators. GJ 3470 b, in contrast, shows significant features in its HST/WFC3 transmission spectrum (Benneke et al. 2019b), indicating that its clouds are lower in the atmosphere, more patchy, or some combination of the two. The reason for the differing levels of cloudiness between these three - otherwise similar - exoplanets is not clear at this time.  \nThe short period of the planetary orbit (3.337 days) places GJ 3470 b at the outer edge of the Neptune Desert, indicating that it may have had a dynamic migration and evolutionary history at earlier times. Transit observations of the Lyα (Bourrier et al. 2018) and  \nHe I (Palle et al. 2020) lines show that the planet is currently losing ∼ 10 10 g s -1 into space. This massloss rate would imply that GJ 3470 has lost roughly 40% of its mass over the 2 Gyr system lifetime (Bourrier et al. 2018). Rossiter-McLaughlin measurements of the spin-orbit alignment of the system (Stef'ansson et al. 2022) show that GJ 3470 b is on a nearly-polar orbit ( λ = 98 +15 -12 deg.). Spitzer secondary eclipse observations (Benneke et al. 2019b) indicate the orbit is not circular, and coupled with radial velocity observations the most recent analyses give e = 0 . 125 ± 0 . 04 (Kosiarek et al. 2019; Stef'ansson et al. 2022).  \nThe combination of the present-day high mass-loss rate, polar orbit, and eccentric orbit point towards the scenario that a more massive version of GJ 3470 b was dynamically scattered inwards towards the star GJ 3470 at some point after the planet formed (Stef'ansson et al. 2022), and subsequently lost a significant fraction of its outer envelope (Palle et al. 2020). The planet as we see it today is likely in the final stages of tidal circularization and atmospheric mass-loss on its way 'out' of the Neptune Desert at the end of this turbulent migration history.  \nMeasurements of GJ 3470 b's present-day atmospheric carbon, oxygen, and sulfur abundances are one method to constrain the planet's initial formation location Turrini et al. (2021); Crossfield (2023), and to better understand the formation and migration histories of subNeptunes.", '2. OBSERVATIONS, DATA ANALYSIS, AND FITTING': "In our analysis of the transmission spectrum of GJ 3470 b, we used a combination of new data collected by JWST/NIRCam and archival transit observations from HST/WFC3 and Spitzer/IRAC. The WFC3 transit observations and some Spitzer data have been previously analyzed and reported on in Benneke et al. (2019b) For this work, we re-extracted and re-fit all of the archival data to minimize possible offsets in the overall transit depths between the different observations. These offsets are often caused by using different priors on the planetary orbital parameters and stellar limb-darkening coefficients, and we wished to fit all of the available data using a single set of orbital parameters and a set of limbdarkening coefficients drawn from a single stellar model.  \nIn addition to the transit observations, we used the Gaia parallax (Gaia Collaboration et al. 2023) for GJ 3470 to fit the stellar spectral energy distribution (SED) and, hence, the stellar radius. Together with the stellar density we measure via the transit lightcurves and previous radial velocity observations (Kosiarek et al. 2019),  \n/uni00000029/uni00000016/uni00000015/uni00000015/uni0000003a/uni00000015  \n/uni00000029/uni00000017/uni00000017/uni00000017/uni0000003a  \nFigure 1. The spectroscopic and broadband NIRCam lightcurves of GJ 3470 b. The top panels show the raw F322W2 and F444W lightcurves binned in the spectral direction into 0.015 µ m-wide spectral bins, with the transit of GJ 3470 b clearly visible at all wavelengths. The bottom panels show the broadband F322W2 and F444W lightcurve data, without their associated error bars for visual clarity.  \n<!-- image -->  \nthis allowed us to make new estimates for the absolute radii and masses for the star GJ 3470 and the planet GJ 3470 b. We used the values we determined for GJ 3470 b's radius and mass in our atmospheric retrievals.  \nWe also update the results from long-term photometric monitoring of the variability of GJ 3470. Previous results estimated the star GJ 3470's rotation period and variability amplitude (Biddle et al. 2014) based on one year of monitoring data. We now have eleven years of photometric monitoring data from the Tennessee State University automated imaging telescope, which allowed us to refine the rotation period measurement and the variability amplitude at this period. As described below, based on previous analyses (Benneke et al. 2019b), the stellar variability we observe does not significantly influence our results.  \nWe briefly summarize all of these steps here. Complete details and a fuller description of the data reduction, lightcurve fitting, and stellar SED and variability analyses are given in Appendices A and B.  \n2.1. Observations and Data Reduction 2.1.1. JWST NIRCam  \nWe observed two transits of GJ 3470 b using JWST/NIRCam as a part of GTO program 1185 (PI Thomas Greene, Greene et al. 2017). We observed the first transit on UT 31 October 2022, using the F322W2 filter together with the grismr disperser on NIRCam's longwave (LW) detector, and the F210M filter together with the wlp8 weak lens array on NIRCam's shortwave (SW) detector. For the F322W2 transit, we configured both detectors to use the subgrism256 subarray in bright2 readout mode with four readout amplifiers and using four groups per integration. Both the LW and SW detectors took 1681 integrations in the complete exposure. In total, we observed the first transit of GJ 3470 b for 5.67 hours.  \nWe observed the second transit of GJ 3470 b on UT 3-4 November 2022 using the F444W filter on the LW detector. The LW and SW detectors were configured similarly to the F322W2 transit, except that for this second transit, we increased the number of groups per integration to seven and decreased the number of integrations per exposure to 1009. This longer integration time kept the peak counts from GJ 3470 in the F444W  \nFigure 2. Three independent reductions of the NIRCam F322W2 and F444W data agree within their uncertainties The top panel shows the F322W2 and F444W transmission spectra, as reduced using the Pegasus , Eureka , and tshirt pipelines, plotted with their associated 1 σ uncertainties using constant 0.015 µ m-wide spectral bins. All of the F322W2 spectra have been offset down by 163 ppm to correct for a bias offset in the F322W2 transit depths compared to the F444W and Spitzer data (see Appendix B.3). All three reductions agree on the shape of the transmission spectrum and show features due to H 2 O (2.5-3.2 µ m), CH 4 (3.2-3.5 µ m), SO 2 (3.9-4.1 µ m), and CO 2 (2.6-2.9 µ m and 4.2-4.5 µ m). The absorption cross sections of the prominent species in this wavelength range, many of which are visually identifiable in the bottom panel.  \n<!-- image -->  \nfilter roughly the same as in the F322W2 data. Again, we observed this second transit for 5.67 hours.  \nOur fiducial NIRCam data reduction used the Pegasus pipeline (https://github.com/TGBeatty/ PegasusProject), which was initially developed for the NIRCam commissioning observations. To ensure the reproducibility of our results, we performed two additional analyses of the NIRCam spectra using the Eureka! (Bell et al. 2022), and tshirt (Ahrer et al. 2023; Bell et al. 2023) data reduction pipelines. Using the same system parameters, limb-darkening coefficients, and a freely-fit linear background trend, we find strong agreement between all three spectra, with differences of χ 2 / dof = 0 . 13 for Pegasus vs. Eureka! and χ 2 / dof = 0 . 21 for Pegasus vs. tshirt . The three spectra all show the same spectroscopic features with no major systematic differences between them (Figure  \n- 2). We chose the Pegasus spectrum as our fiducial case for atmospheric modeling due to its ability to reduce the archival HST and Spitzer data. The complete details of all three reduction pipelines and their results are described in Appendices A and B.", '2.2. HST WFC3': 'HST observed three transits of GJ 3470 b using WideField Camera 3 (WFC3) and the G141 disperser as a part of GO program 13665 (PI Benneke) using four orbits of HST time per transit (Benneke et al. 2019b). To provide consistency with our analysis and assumed system and orbital parameters, we re-reduced and re-fit two of the transit observations: those taken on UT 13 March 2015 and UT 22 October 2015. We reduced the spectroscopic observations and extracted broadband and spectroscopic lightcurves using previously described tech-  \nniques Beatty et al. (2017). The broadband lightcurve covered the wavelength range from 1.12 µ m to 1.66 µ m in a single channel, while the spectroscopic lightcurves covered this same wavelength range in fifteen bins that were each 0.036 µ m wide. This wavelength range and the spacing match the analysis of these data in (Benneke et al. 2019b).', '2.3. Spitzer IRAC': 'The GJ 3470 system has been repeatedly observed in transit and eclipse by Spitzer at 3.6 µ m and 4.5 µ m. The combined results from three of the 3.6 µ m transits and three of the 4.5 µ mtransits have also been described in Benneke et al. (2019b). Here, we reduced and simultaneously fit a total of seven 3.6 µ m and eight 4.5 µ m Spitzer transit observations of GJ 3470 (which includes the three in each channel that are already in the literature). These data were taken as a part of Spitzer GO programs 90092 (PI Desert), 80261 (PI Demory), and 13140 (PI Kreidberg). All of these observations used the same instrument and exposure time settings in subarray mode with 0.4-second exposures and PCRS peak-up. We reduced these data and extracted photometry using previously described procedures (Beatty et al. 2018, 2019).', '2.4. Broadband and Spectroscopic Lightcurve Fitting': "We first performed a simultaneous fit to broadband HST, broadband JWST, and Spitzer transit lightcurves to refine the general orbital and planetary parameters of the GJ 3470 system. Our results from this joint fit to the broadband photometry (Table 2, Figure 6) were generally consistent with previous measurements (Biddle et al. 2014; Awiphan et al. 2016; Benneke et al. 2019b; Kosiarek et al. 2019). We did find that the planetary orbit was slightly less eccentric and the transit depth was mildly shallower (Benneke et al. 2019b; Kosiarek et al. 2019), due to our analysis with the new NIRCam data measuring GJ 3470 b's impact parameter to be b = 0 . 31 ± 0 . 02, which is slightly smaller than those older estimates of b ≈ 0 . 4 (Dragomir et al. 2015; Awiphan et al. 2016; Benneke et al. 2019b; Kosiarek et al. 2019).  \nWe measured the spectroscopic transit depths of GJ 3470 b using the JWST/NIRCam observations and the archival HST/WFC3 data. For these fits, we fixed all of the transit properties except for the planet-to-star radius ratio to the values determined from the joint photometric fit. The JWST data was relatively clean and preferred a simple linear background temporal trend to remove systematic trends in the lightcurves. The HST/WFC3 data showed the usual hook-like systematic  \ntrends. We fit these by trending the data within each orbit using a hooked exponential and a visit-long linear slope, the details of which are described in more detail in Beatty et al. (2017). We found the mean depth of the WFC3 spectrum to be mildly shallower than previously reported in Benneke et al. (2019b), a result of more refined measurements of a/R ∗ and cos( i ) driven by the JWST data. As described in Appendix B.3, we found that the F322W2 spectrum was offset 163 ± 34 ppm high relative to the F444W and Spitzer transit depths. Welbanks et al. (2024) recently found a roughly similar offset between F322W2 and F444W transmission spectra of WASP-107 b, and Xue et al. (2024) marginally detected an offset in their F322W2 and F444W transmission spectra of HD 209458 b. Possibly, the offset is due to flux-dependent non-linearity effects in the F322W2 observations in all three data sets - though more investigation is necessary to confirm the cause of the offsets between F322W2 and F444W transit data.", '2.5. Stellar SED Modeling': "We used a set of catalog 2MASS JHK (Skrutskie et al. 2006) and AllWISE W1 and W2 (Cutri & et al. 2014) magnitudes for the GJ 3470 system to fit a spectral energy distribution (SED) model to the stellar emission. In conjunction with the Gaia (Gaia Collaboration et al. 2023) parallax for the system, this allowed us to estimate a stellar radius for the star GJ 3470. To model the SED, we used BHAC15 spectra (Baraffe et al. 2015) for the stellar SED.  \nThis stellar SED fitting (Figure 8) allowed us to measure the radius of the star GJ 3470 to be R ∗ = 0 . 50 ± 0 . 01 R ⊙ . We then used the mean value of R P /R ∗ as measured in our joint broadband transit fits to estimate the planetary radius of GJ 3470 b to be R p = 4 . 22 ± 0 . 09 R ⊕ . Using the stellar density we measure as a part of the joint broadband transit fits, we constrain GJ 3470's stellar mass to be M ∗ = 0 . 44 ± 0 . 04 M ⊙ , and using a mass-ratio measured via radial velocity observations (Kosiarek et al. 2019), we estimate a planetary mass of M p = 11 . 1 ± 0 . 9 M ⊕ . Our atmospheric retrievals use these absolute stellar and planetary radii and masses.", '2.6. Stellar Variability Monitoring': "Previous photometric monitoring of GJ 3470 during the 2012-2013 season showed that the star itself is slightly variable at the presumed stellar rotation period of 21 days (Biddle et al. 2014). This same monitoring campaign now spans from 2012 to 2023 and has given us a more precise measurement of stellar variability (Figure 9).  \nUsing these new additional observations, we find that the GJ 3470 system shows sinusoidal brightness modu-  \nFigure 3. The combined NIRCam and HST/WFC3 transmission spectrum of GJ 3470 b showing major absorption features. The measured transit depths are plotted with their associated 1 σ uncertainties. The NIRCam data have been binned into 0.06 µ m-wide spectral channels (black points) for visual clarity. The background grey points show the NIRCam data binned in 0.015 µ m-wide channels, which is the direct result of our analysis and was used in our atmospheric modeling. The solid dark-grey line shows the best fit model transmission spectrum from the atmospheric free retrieval analysis. The colored regions illustrate the contributions of various molecules and clouds to the best fit model. Note that the new NIRCam data and the refined HST transit depths strongly disfavor the previously suggested Mie scattering feature near 3 µ m (dashed grey line, Benneke et al. 2019b).  \n<!-- image -->  \nlations with an amplitude of 2545 ± 130 ppm at a period of P = 21 . 6239 ± 0 . 0014 days. This variability amplitude is roughly half the value measured previously using only the 2012-2013 data (Biddle et al. 2014) and roughly half the variability amplitude used in previous modeling of how stellar variability could affect GJ 3470 b's transmission spectrum (Benneke et al. 2019b).  \nPrevious work (Benneke et al. 2019b) modeled the wavelength dependence of GJ 3470's variability and its effect on transmission spectra using a variability amplitude approximately twice that which we measured and found that it was at most 20 ppm in the WFC3 bandpass and < 10 ppm over the F322W2 and F444W wavelength range. Since even our conservative estimate for the effect of stellar variability - assuming a wavelengthindependent amplitude - is on the order of (or smaller) than our estimated depth uncertainties and the interinstrument offsets, we do not expect the variability seen in the GJ 3470 system to affect our results significantly.", '3. ATMOSPHERIC MODELING': "We interpreted the transmission spectra of GJ 3470 b using one-dimensional atmospheric models spanning a wide range of physico-chemical assumptions. First, we  \nuse parametric models coupled with Bayesian parameter estimations using nested sampling (Buchner et al. 2014) (e.g., 'free' atmospheric retrievals, see Methods) to quantify the detection of the multiple molecular features present in the spectrum of the planet. This paradigm independently fits for the chemical abundances, vertical pressure-temperature structure, and clouds in the planetary atmosphere without any assumptions of radiativeconvective thermo-chemical equilibrium.  \nThe Bayesian model comparison results in strong detections of H 2 O (6.3 σ ),CO 2 (7.3 σ ), CH 4 (3.8 σ ), and SO 2 (4.0 σ ), confirming the importance of the spectral features visible in the planetary spectrum (Figure 3). Although the retrieved molecular abundances can be sensitive to the presence of possible offsets between the different instruments in GJ 3470 b's spectrum, the reported detections are robust regardless of any instrumental shifts and do not depend on the 163 ppm offset between NIRCAM F322W2 and F444W. The retrieved atmospheric abundances are suggestive of a super-solar metallicity with estimates of ≳ 100 × solar (e.g., log 10 ( X H 2 O ) = -1 . 08 +0 . 43 -0 . 52 , log 10 ( X CH 4 ) = -4 . 05 +0 . 25 -0 . 27 , log 10 ( X CO 2 ) = -2 . 47 +0 . 61 -0 . 43 , and log 10 ( X SO 2 ) = -3 . 57 +0 . 26 -0 . 25 ).  \nFigure 4. A comparison of the results from the 1D-RCPE and the free-retrieval atmospheric modeling, demonstrating their agreement. The free-retrieval's posterior distributions (gray) are shown against the abundance profiles (colored lines) from the 1D-RCPE grid modeling. The dotted lines show the expected composition of a Solar metallicity and Solar C/O atmosphere. The dash-dotted lines show the expected abundance profiles for a 125 × solar metallicity atmosphere with a C/O of 0.35, as estimated for GJ 3470 b, under chemical equilibrium. The solid colored lines and shaded regions show the median 1 σ and 2 σ profiles for GJ 3470 b under the disequilibrium conditions in the 1D-RCPE modeling. Note that for all molecules except for CH 4 and SO 2 the abundances estimated from the 1D and free-retrieval modeling match the equilibrium abundances of the 125 × solar metallicity atmosphere with a C/O of 0.35. The indicates that CH 4 is likely being depleted and SO 2 enhanced via disequilibrium processes.  \n<!-- image -->  \nThen, to obtain a self-consistent estimate of the atmospheric metallicity and carbon-to-oxygen ratio we use 1-dimensional radiative-convective-photochemical equilibrium models (1D-RCPE, Bell et al. 2023), which we describe in more detail in Appendix E. Briefly, the 1DRCPEmodels self-consistently compute the atmospheric pressure-temperature profiles and gas volume mixing profiles using given values for the incident stellar flux, internal temperature, elemental abundances (via a scaled metallicity and C/O), and the vertical mixing strength. Using the same Bayesian inference methodology we derive the temperature, gas, and molecular weight profiles derived from the pre-computed 1D-RCPE models above are interpolated and post-processed via a transmission geometry radiative transfer routine (see Methods). As  \nwith the 'free' retrieval the atmospheric models account for the effects of clouds and hazes using parameterizations for their inhomogeneous limb cover (see Methods).  \nThe 1D-RCPE models agree with the free-retrieval interpretation of a super-solar atmospheric metallicity for GJ 3470 b, retrieving a metallicity of 125 ± 40 times Solar (M / H = 2 . 1 ± 0 . 15 dex, Lodders et al. 2009), and a slightly sub-Solar carbon-to-oxygen ratio of C / O = 0 . 35 ± 0 . 1. The inferred molecular profiles (e.g., molecular volume mixing ratios as a function of height) are generally consistent with the inferred volume mixing ratios from the constant with height free-retrieval as shown in Figure 4. The general agreement in the atmospheric composition between the flexible parametric free-retrieval and the 1D-RCPE retrieval confirms the  \n<!-- image -->  \nFigure 5. The impact of composition and temperature on the SO 2 abundance in cool ( < 800 K) sub-Neptune atmospheres. The left panel illustrates the column-averaged (1 µ bar - 0.1 mbar) volume mixing ratio of SO 2 (contours) as a function of the atmospheric metallicity (relative to Solar) and the carbon-to-oxygen ratio for the nominal temperature of GJ 3470 b. There is a steep fall off in SO 2 with increasing C/O and increasing carbon-to-oxygen ratio. The right panel shows the dependency of the column-averaged SO 2 abundances with planetary irradiation temperature for representative metal enrichment. The presence of SO 2 at these temperatures substantially restricts the possible atmospheric composition. The retrieved SO 2 abundance (navy contour on the left and the shaded box on the right) strongly indicates an atmosphere with ∼ 100 × stellar metallicity and a sub-stellar carbon-to-oxygen ratio.  \n<!-- image -->  \ninterpretation of a super-solar metallicity atmosphere affected by disequilibrium processes (e.g., photochemistry) that give rise to SO 2 production. As shown in Figure 3, this atmosphere model adequately explains the observed panchromatic spectrum, capturing nearly all the salient features.", '4. DISCUSSION': "The detection of SO 2 in the atmosphere of GJ 3470 b extends the regime for sulfur-based photochemistry to lower planetary masses and cooler temperatures than previously probed. SO 2 has been detected in the atmospheres of WASP-39 b (90 M ⊕ , T eq = 1166K, Alderson et al. 2023; Rustamkulov et al. 2023) and WASP107 b (30 M ⊕ , T eq = 770K, Dyrek et al. 2024). GJ 3470 b is significantly lower in mass (11 M ⊕ ) and at a significantly lower equilibrium temperature (T eq = 600 KBonfils et al. (2012)) than either of these planets. The detection of SO 2 in GJ 3470 b's atmosphere opens the door to using sulfur abundances to recover the formation history of sub-Neptunes (Turrini et al. 2021; Crossfield 2023).  \nCurrent models of sulfur photochemistry (Tsai et al. 2021; Polman et al. 2023; Tsai et al. 2023) indicate that at GJ 3470 b's temperature and mass, the SO 2 abundance should be low and more difficult to detect. This is predicted to occur because the oxidation reactions necessary to create SO 2 from H 2 S and OH would proceed over an order-of-magnitude more slowly than in the atmosphere of WASP-39b at the same metallicity and C/O due to the lower temperature on GJ 3470 b (Zahnle et al.  \nTable 1. Free-retrieval Atmospheric Properties.  \n2016; Tsai et al. 2023). Note, however, that the existing literature on sulfur photochemistry predominantly focuses on Jupiter-like planetary atmospheres with metallicities of ∼ 1-10 × Solar.  \nThe presence of detectable SO 2 in GJ 3470 b's atmosphere likely indicates that the planetary atmosphere  \nhas a much higher metallicity, resulting in abundant water vapor. Models of sulfur photochemistry in Tsai et al. (2023) predict that if the planetary atmosphere has a large amount of water vapor - consistent with the high H 2 O abundance we measure - this will allow the UVdriven photolysis of H 2 O on the planetary dayside to generate large amounts of OH radicals and overcome the slower chemical reaction rates due to the lower temperature. As shown in Figure 5, our modeling shows that to generate the SO 2 abundance we observe in GJ 3470 b (log 10 ( X SO 2 ) = -3 . 55 +0 . 30 -0 . 52 ) the atmosphere needs to have a metallicity of ∼ 10-300 × Solar and a C/O ratio ≲ 0 . 35, consistent with other parts of our analysis.  \nThe H 2 Oand CO 2 abundances we measure in GJ 3470 b's atmosphere contrast with the low inferred CH 4 abundance of log 10 ( X CH 4 ) = -4 . 05 +0 . 25 -0 . 27 . This is significantly higher than the 3 σ upper-limit on the methane abundance of log 10 ( X CH 4 ) < -4 . 9 reported in Benneke et al. (2019b). At the same time, our measured CH 4 abundance is approximately 10 × to 100 × lower than one would expect from thermochemical equilibrium calculations at the inferred metallicity and C/O of the atmosphere (Figure 4). UV photolysis is unlikely to cause the relatively low CH 4 depletion at the pressures probed by these data (Miller-Ricci Kempton et al. 2012; Moses 2014). It is possible that the CH 4 level in GJ 3470 b's atmosphere is being reduced by vertical chemical quenching and the transport of CH 4 -poor and CO-rich gas up to the pressures probed by these observations Venot et al. (2014). CH 4 depletion via vertical disequilibrium processes has been suggested for other exoplanet atmospheres recently observed with JWST (Rustamkulov et al. 2023; Xue et al. 2024). In GJ 3470b, the CH 4 abundance may be further reduced due to tidal heating. Spitzer eclipse observations show that GJ 3470 b's orbital eccentricity is significantly non-zero, with e cos ω = 0 . 01457 ± 0 . 0075 (Benneke et al. 2019b), consistent with the RV observations (e.g., Kosiarek et al. 2019; Stef'ansson et al. 2022). For a Neptune-like tidal quality factor of Q = 10 4 , the planet's internal temperature would be T int ≈ 550 K. This would increase vertical mixing and decrease the CH 4 abundance to match our observations (Fortney et al. 2020). Note that our modeling only provides a loose constraint of T int = 250 ± 200 K.  \nThe strong signature of SO 2 in our data supports the interpretation that GJ 3470 b possesses a high metallicity, near Solar C/O, atmosphere with reduced CH 4 from disequilibrium processes. In principle, an analysis using atmosphere models assuming thermochemical equilibrium might find a consistent solution that accounts for the H 2 O, CO 2 , and CH 4 abundances we measure - likely  \na Solar metallicity and very low C/O atmosphere as previously reported for GJ 3470 b in Benneke et al. (2019b) However, the significant SO 2 we measure in GJ 3470 b is inconsistent with equilibrium models (Tsai et al. 2021; Polman et al. 2023). As discussed above, the presence of SO 2 is itself indicative of disequilibrium processes and a high metallicity atmosphere (see also Figure 5). Thus, we consider this the preferred scenario.  \nRecently, Turrini et al. (2021) and Crossfield (2023) have proposed using the combination of carbon, oxygen, and sulfur abundances in exoplanet atmospheres to estimate their initial formation locations in the protoplanetary disk, and possibly the formation process itself (i.e. planetesimal vs. pebble accretion, Crossfield 2023). Specifically, both papers suggest that the C / S and O / S ratios in a planetary atmosphere should increase by a factor of 3x-20x relative to Solar as the initial formation location increases from 5-30 AU in the case of pebble accretion. The atmospheric abundances that we measure in GJ 3470 b yield rough estimates of C / S = 15 ± 10 times Solar and O / S = 20 ± 10 times Solar. This suggests pebble accretion at 10-30 AU from the host star, but better constraints are needed on the CO abundance and the H 2 S abundance to say anything definite.  \nThe results from these observations further illustrate the new complexity in exoplanet atmospheres that is becoming visible with JWST observations. The observed CH 4 depletion and strong signature of SO 2 indicate that GJ 3470 b's atmosphere is a dynamic system - both kinetically and chemically - that is not easily summarized without an accounting of all the major carbon- and oxygen-bearing molecules. Though the high metallicity and Solar C/O we infer for GJ 3470 b are both more in line with the theoretical expectations from planet formation models, the compositional diversity being revealed by JWST data indicates that relating observations of atmospheric compositions to formation histories may be a more idiosyncratic process than has been expected. The clear detection of SO 2 on a cool planet like GJ 3470 b also indicates that photochemistry is an important consideration across a wide range of temperatures. Scheduled future observations of GJ 3470 b's transmission spectrum using JWST/NIRSpec G395H will hopefully further illustrate the intriguing complexity of this atmosphere and further constrain its formation history by providing an independent measurement of the atmospheric SO 2 abundance.  \nThis work has made use of NASA's Astrophysics Data System, the Extrasolar Planet Encyclopedia at exoplanet.eu (Schneider et al. 2011), the SIMBAD database operated at CDS, Strasbourg, France (Wenger et al. 2000), and the VizieR catalog access tool, CDS, Strasbourg, France (Ochsenbein et al. 2000). This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. L.W. acknowledges support for this work provided by NASA through the NASA Hubble Fellowship grant HST-HF2-51496.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS526555. T.J.B. and T.P.G. acknowledge funding support from the NASA Next Generation Space Telescope Flight Investigations program (now JWST) via WBS 411672.07.05.05.03.02. M.R.L. acknowledges NASA XRP award 80NSSC19K0446 and STScI grant HSTAR-16139. M.R.L. and L.W. acknowledge Research Computing at Arizona State University for providing HPC and storage resources that have significantly contributed to the research results reported within this manuscript. K.O. acknowledges support from JSPS KAKENHIGrant Number JP23K19072. M.M. acknowledges funding from NASA Goddard Spaceflight Center via NASA contract NAS5-02105. G.W.H. acknowledges long-term support for our automatic telescopes from NASA, NSF, and Tennessee State University. We thank Marcia Rieke for allocating the NIRCam time for this program.  \nFacilities: HST (WFC3), JWST (NIRCam), Spitzer (IRAC)", 'A.1. JWST Data Reduction': 'As mentioned, to ensure the reproducibility of our results we performed three separate analyses of the NIRCam spectra using the Pegasus , Eureka! (Bell et al. 2022), and tshirt (Ahrer et al. 2023; Bell et al. 2023) data reduction pipelines (Figure 2). We reduced both of the LW spectroscopic transit observations using three separate data reduction pipelines, which we then compared as a self-consistency check. We performed a single reduction on the SW photometric observations.', 'A.1.1. Pegasus Data Reduction': 'Our fiducial NIRCam data reduction used the Pegasus pipeline (https://github.com/TGBeatty/ PegasusProject), which was initially developed for the NIRCam commissioning observations. Our reduction began with background subtraction on the rateints files provided by version 1.10.2 of the jwst pipeline using CRDS version 11.17.0. For each rateints file, we first performed a general background subtraction step by fitting a two-dimensional, second-order spline to each integration using the entire 256 × 2048 rateints images. To fit the background spline, we first masked out image rows 5 to 75 to not self-subtract light from the GJ 3470 system and from rows 130 to 140 in the F322W2 images to mask a nearby background star. We then performed a single round of 3 σ clipping on the unmasked portions of the image. We then individually fit a spline to each amplifier region on the image using the unmasked, unclipped pixel values with a median-box size of 20 pixels. Fitting the background spline on a per-amplifier basis is necessary to remove residual bias offsets between the amplifier regions. We extrapolated the combined background spline for the whole image over the masked portions near GJ 3470 and subtracted it from the original image values. Visual inspection of the rateints images showed that in roughly 5% of the integrations the reference pixel correction failed for at least one of the amplifier regions, so after the spline fitting and subtraction we re-ran the reference pixel correction using hxrg-ref-pixel (https://github.com/JarronL/hxrg ref pixels). Finally, we attempted to remove some of the red-noise caused by the NIRCam readout electronics along detector rows, by calculating the biweight mean of each row using pixels from column 1800 onwards for the F322W2 images and up to column 600 for F444W (both chosen to avoid light  \nfrom GJ 3470), and then subtracting this mean from each row. Visually, this removed most of the horizontal banding typical for NIRCam grismr images.  \nWe then extracted broadband and spectroscopic lightcurves from our background-subtracted images. To do so, we fit the spectral trace using a fourth-order polynomial and then used optimal extraction to measure the 1D spectrum in each image. We performed three rounds of iterative profile estimation for the optimal extraction routine, after which we judged the profile fit to have converged. Using the resultant 1D spectra, we extracted a broadband lightcurve from 2.45 µ m to 3.95 µ m at F322W2 and from 3.89 µ m to 4.97 µ m at F444W. For the spectroscopic lightcurves, we subdivided each of these wavelength regions into individual 0.015 µ m-wide spectral channels. For both the broadband and spectroscopic lightcurves, we linearly interpolated over each spectral column to account for partial-pixel effects in the wavelength solutions.', 'A.1.2. Eureka! Data Reduction': 'We also reduced the observations using the opensource Eureka! pipeline (Bell et al. 2022). Our Eureka! reduction of the F322W2 data followed the same method as described by Bell et al. (2023), with the only difference being the package versions used. Our Eureka! reduction used version 0.10.dev15+g80126b56.d20230613 of the Eureka! pipeline, CRDS version 11.17.0 and context 1093, and jwst package version 1.10.2. The Eureka! Control Files and Parameter Files we used are available for download (https://doi.org/10.5281/ zenodo.10779291). In summary, we cropped the full subarray images only to use y-pixels 5-64 and x-pixels 151709 and then corrected for the curvature of the spectral trace. We then performed per-column, per-integration background subtraction on the Stage 2 outputs; and performed optimal spectral extraction using the pixels within 9 pixels of the center of the flattened spectral trace. For each spectroscopic lightcurve, we then masked 10-sigma outliers compared to a boxcar-filtered version of the lightcurve (with a boxcar width of 50 integrations) to remove any remaining cosmic rays while not removing sharp features like transit ingress/egress.  \nOur Eureka! reduction of the F444W data used the same package versions as the F322W2 observations. We ran Stages 1 and 2 of Eureka! using the same parameters as for the F322W2 reduction. For Stage 3, we cropped each image to only include x-pixels 738-2000 (since the spectra fall on a different region of the de-  \ntector), and we used a smaller aperture half-width of 4 when extracting the spectra.', 'A.1.3. tshirt Data Reduction': 'Finally, we also similarly reduced the data as for WASP-80 b (Bell et al. 2023) using the tshirt (Ahrer et al. 2023; Bell et al. 2023) pipeline. We used data from SDP VER = 2022 3b, JWST version 1.6.0, CRDS version 11.16.5 for the stage 1 ( calwebb detector1 ) processing, with some modifications. The CRDS context was jwst 1011.pmap for the F322W2 filter and jwst 1009.pmap for the F444W filter, but the reference files used were the same. The modifications to the calwebb detector1 were to replace the reference pixel step with a Row-by-row Odd/Even Subtraction (ROEBA) step using sky pixelsSchlawin et al. (2023). We also skip the dark current step as the subarray dark subtraction adds noise to the time series, and finally we adjust the jump step rejection threshold to be 6 σ to avoid spurious false positive jump step detections. At the completion of calwebb detector1 , we divide the NRCALONG detector image by the imaging flat field in the corresponding band.  \nFor the spectroscopic extraction, we subtracted the background column-by-column with a linear robust fit and summed over a rectangular aperture with a covariance-weighted extraction profile (Schlawin et al. 2020). The background linear fit was from pixels that were 10 to 30 pixels from the source centroid. The aperture was a 10 pixel wide aperture centered on pixel Y=34 for the F322W2 filter and Y=31 for the F444W filter. We assumed a read noise correlation of 0.08 between pixels and a read noise of 14 e -.', 'A.1.4. Shortwave F210M NIRCam Photometry': "NIRCam can collect shortwave photometric observations simultaneously with longwave spectroscopic observations. To reduce both of the photometric F210M observations taken along with the spectroscopic F322W2 and F444W transit observations, we used the Eureka! data analysis pipeline version 0.9, (Bell et al. 2022) with CRDS version 11.16.20, jwst version 1.6.2. The F322W2 reduction used the CRDS context 1041 and the F444W reduction used 1046. We began by running Stages 1 to 4 of Eureka! on the uncalibrated, uncal , image files. For Stage 1, we used the default jwst settings, except that we changed the cosmic ray jump rejection threshold to 10 σ to avoid excessive false positive cosmic ray rejections. We ran Stage 2 with default settings, except that we did not perform the photometric calibration step.  \nTo run Stage 3 on the F210M data, we first cropped the images to focus on the source location on the de-  \ntector. For the photometry collected along with the F322W2 transit, we cropped to y-pixels 0-255 and xpixels 550-1450, and y-pixels 20-250 and x-pixels 10002040 for the F444W transit. We then used a cubic interpolation method to fill in bad pixels and outliers in these cropped images. We identified bad pixels using the data quality array, and 7 σ outlier pixels along the time axis. We then used Eureka! 's version of the ROEBA (Schlawin et al. 2023) background correction technique to remove 1/ f noise in the images. To do so, we excluded pixels within 350 pixels of the location of GJ 3470 on the detector. We also measured the centroid of GJ 3470's point-spread function in each integration in a two-step process. We first estimated the centerof-light over an entire cropped image to get an initial guess at the centroid location. We then performed another center-of-light estimate using just an 11-pixel by 11-pixel box centered on the initial centroid guess.  \nWe extracted lightcurves using the resultant calibrated images from both F210M photometric observations using simple aperture photometry. We used a circular extraction aperture with a radius of 65 pixels for both observations, and a background annulus with inner and outer radii of 70 and 90 pixels for F322W2's shortwave observations and 75 and 90 pixels for F444W's shortwave observations. Finally, in Eureka! 's Stage 4, we clipped 8 σ outliers compared to a boxcar smoothed version of the data using a boxcar width of 50 integrations. The resulting F210M lightcurves contained 1681 points for the F322W2 transit with a median uncertainty of ± 234 ppm per point, and 1009 points for the F444W transit with a median uncertainty of ± 170 ppm per point.", 'A.2. HST Data Reduction': "We reduced the spectroscopic observations and extracted broadband and spectroscopic lightcurves using the pipeline described in Beatty et al. (2017). In brief, we extracted 1D spectra from each of the scanned, upthe-ramp readouts using the G141 disperser after performing a set of bad-pixel and scan-angle corrections. We summed the 1D spectra from each of these subexposures together to make a single 1D spectrum for each exposure. Next, we extracted a broadband lightcurve and a set of spectroscopic lightcurves using the summed 1D spectra. The broadband lightcurve covered the wavelength range from 1.12 µ m to 1.66 µ m in a single channel, while the spectroscopic lightcurves covered this same wavelength range in fifteen bins that were each 0.036 µ m wide. This wavelength range and the spacing match the previous analysis of these data (Benneke et al. 2019b). The broadband and spectroscopic HST/WFC3  \nlight curves for each visit contained 72 points split evenly between the visit's four HST orbits. The median flux uncertainty for the broadband lightcurve was 100 ppm, and the median flux uncertainty among all the spectroscopic lightcurves was 286 ppm.", 'A.3. Spitzer Data Reduction': "We began our data reduction and photometric extraction process from the basic calibrated data (BCD) images. We first subtracted the background from the images and corrected for bad pixels. We then used the background-subtracted, bad-pixel corrected images to measure the pixel position of GJ 3470 in each image using a two-dimensional Gaussian. Note that we used these corrected images only to estimate the background and to measure the position of GJ 3470 - we used uncorrected background-subtracted images for the photometric extraction.  \nWe extracted raw photometry for GJ 3470 using a circular extraction aperture centered on the star's position in each image. We used an aperture radius of 2.4 pixels. For reference, the average full-width half-maximum of GJ 3470's point spread function was 2.05 pixels at 3.6 µ m and 2.00 pixels at 4.5 µ m. Since the fitting process for these observations was time-intensive, we did not perform a complete optimization to determine the best extraction aperture size. Instead, we used a fixed aperture radius of 2.4 pixels for all the datasets, matching the average optimum aperture size estimate previously in Benneke et al. (2019b). We performed a limited test of our aperture size by extracting and fitting photometry for aperture radii of 2.2 and 2.6 pixels. In both cases, the log-likelihoods of the resulting best fits were lower, the scatter in the residuals higher, and the transit depths consistent with our optimum aperture at 2.4 pixels.  \nWe trimmed the first 250 data points for each observation to remove residual ramp effects and performed a single round of 5 σ clipping to remove outliers. The median uncertainty at 3.6 µ m was 4462 ppm for all the observations, and at 4.5 µ m was 6073 ppm.", 'B. LIGHTCURVE FITTING': "To measure the transmission spectrum of GJ 3470 b, we first performed a simultaneous, joint fit of all six of the broadband lightcurves (Figure 6) and then used the resulting orbital properties for GJ 3470 b as fixed input to our spectroscopic fits to the WFC3, F322W2, and F444W data. Our intention in doing so, and in refitting transit depths that have been previously reported, was to construct a self-consistent transmission spectrum for GJ 3470 b across all of our data without possible offsets caused by different priors on the planetary orbital  \nparameters and different limb-darkening properties. In particular, we measure the transit depth in the WFC3 data to be approximately 100 ppm shallower than in Benneke et al. (2019b) solely owing to more refined measurements of a/R ∗ and cos( i ) driven by the JWST data.  \nFor the joint broadband fit, we used Gaussian priors on GJ 3470 b's a/R ∗ (Awiphan et al. 2016), orbital inclination (Awiphan et al. 2016), e cos ω (Benneke et al. 2019b), e sin ω (Kosiarek et al. 2019), transit center time (Awiphan et al. 2016), and period (Awiphan et al. 2016). We did not impose priors on the planet-to-star radius ratio in any of the observed bands.  \nIn all of our broadband and spectroscopic fitting, we used limb-darkening coefficients calculated from the Kurucz ATLAS9 odfnew stellar models (Castelli & Kurucz 2003). We did so using the nominal spectroscopic stellar properties for GJ 3470 (Kosiarek et al. 2019).", 'B.1. Joint Broadband Transit Fitting': 'To measure the orbital and system properties of GJ 3470 b, we first fit broadband HST, broadband JWST, and Spitzer transit lightcurves simultaneously (Figure 6). We fit a single set of orbital and system parameters for all seven lightcurves and individual values for R p /R ∗ for each data set and filter. The exceptions to this were the two F210M lightcurves, which we fit using a single shared value of R p /R ∗ for both. We used BATMAN (Kreidberg 2015) transit models and performed an MCMC fit to the broadband data using the emcee (Foreman-Mackey et al. 2013) Python package. We imposed Gaussian priors on some of the transit parameters as described above. We judged the MCMC fitting to have converged ny checking the Gelman-Rubin statistic for each parameter was below 1.1, and by visual inspection of the posterior corner plot. We also checked the Gaussianity of the residuals for each of the broadband lightcurve fits using an Anderson-Darling test, and we did not find any significant non-Gaussianity.  \nWe detrended all three of the JWST broadband datasets (F210M, F322W2, and F444W) using linear trends in time and a free normalization factor. We allowed for different slopes between each of the F210M visits. For the HST broadband lightcurve we detrended the usual systematics seen in WFC3 lightcurves using a combined linear ramp and exponential ramp - with exactly the same functional form as we describe below in the fitting process for the WFC3 spectroscopic lightcurves.  \nThe broadband Spitzer lightcurves showed positiondependent flux variations caused by intrapixel effects that are typical for these observations. We used a linear ramp coupled with a Bi-Linear Interpolated Subpixel Sensitivity (BLISS) map (Stevenson et al. 2012) to de-  \nFigure 6. The broadband lightcurve data and bestfit transit models used in this analysis. To refine the orbital parameters of GJ 3470 b, we simultaneously fit archival broadband HST/WFC3 and Spitzer/IRAC data (top panels) together with the new JWST photometric (F210M) and spectroscopic (F322W2 and F444W) data (bottom panels). This provides significantly more precise measurements of the orbital ephemerides and the orbital separation. We use these results to fix the orbital properties of GJ 3470 b in our spectroscopic fitting. The full results of the joint broadband fit are given in Table 2. Note that the F210M panel shows both F210M lightcurves overplotted on top of each other: one from each of the F322W2 and F444W observations.  \n<!-- image -->  \n/uni00000032/uni00000055/uni00000045/uni0000004c/uni00000057/uni00000044/uni0000004f/uni00000003/uni00000033/uni0000004b/uni00000044/uni00000056/uni00000048  \ntrend the Spitzer lightcurves. Recall that there are a total of seven 3.6 µ m transits and eight 4.5 µ m transits of GJ 3470 b that have been observed by Spitzer. To fit all these data simultaneously, we constructed a single, combined, 3.6 µ m BLISS map from all seven of the 3.6 µ m transits, and a single, combined, 4.5 µ m BLISS map from all eight of the 4.5 µ mtransits. This technique has been used previously (Beatty et al. 2019), since the underlying intrapixel sensitivity variations in both irac channels appear to be stable over multi-year timescales (Beatty et al. 2019; Murphy et al. 2023). We did allow for individual background linear slopes between all fifteen Spitzer transits.  \nThe results from the joint broadband fit are listed in Extended Data Table 1. Note that, although we fit for broadband transit depths using the WFC3, F322W2, and F444W data, we did not use these broadband depths in our subsequent atmospheric modeling.', 'B.2. JWST Spectroscopic Transit Fitting': 'Next, we performed spectroscopic transit to fits the F322W2 and F444W data. To do so, we fixed the transit center time, orbital period, scaled semi-major axis, orbital period, eccentricity, and the argument of periastron to the bestfit values from the joint broadband fit (Table 2). As with our analysis of the broadband data, we used quadratic limb-darkening coefficients estimated from the Kurucz ATLAS9 odfnew stellar models (Castelli & Kurucz 2003). We calculated limb-darkening coefficients for each of the spectral channels, and fixed the limb-darkening coefficients in each channel to these model values using the stellar spectroscopic properties from Awiphan et al. (2016). We experimented with using a non-linear limb-darkening law and allowing the limb-darkening coefficients to float as freely fit parameters and found that this did not significantly affect our results.  \nSimilar to our data reduction and extraction approach, we used three separate analysis pipelines to fit the F322W2 and F444W transmission spectra as a  \nself-consistency check and a single analysis pipeline on the WFC3 data. Note that each of the three analysis pipelines we used for the JWST data fit lightcurves extracted from their corresponding data reduction routine described above. We tested the different transit fitting routines on the lightcurves produced by the other pipelines (e.g., Pegasus lightcurve fitting on Eureka! and tshirt timeseries data) and found no significant differences in the end result.  \nWe experimented using linear, quadratic, and no temporal detrending terms in the F322W2 and F444W spectroscopic fits. We found that a linear temporal trend was strongly preferred by a Bayesian Information Criteria (BIC) analysis over both the quadratic (∆BIC ≥ 6) and no detrend (∆BIC ≥ 8) in all spectral channels in both filters. We, therefore, used a linear trend in time in all three of our transit fitting analyses.  \nThe measured transmission spectra from our fiducial Pegasus reduction of the F322W2 and F444W data are given in Tables 3 and 4. Note that in Table 3, the listed depths have not been offset by the 163 ppm offset we estimated for the F322W2 data and discuss in Section 2.4 and Appendix B.3.', 'B.2.1. Pegasus': "We used a BATMAN (Kreidberg 2015) transit model in our Pegasus lightcurve fitting. Most of the transit model parameters were fixed to the values determined by the joint broadband fit, which left the free parameters in our spectroscopic lightcurve fitting to be the planet-tostar radius ratio and the slope and normalization of the background linear trend. We did not impose a prior on any of these parameters. We fit each spectral channel individually.  \nTo fit the spectroscopic lightcurves we performed an initial likelihood maximization using a Nelder-Mead sampler followed by MCMC likelihood sampling. We used the maximum likelihood point identified by the Nelder-Mead maximization as the starting locus for initializing the MCMC chains. To perform the MCMC runs, we used the emcee (Foreman-Mackey et al. 2013) Python package using 12 walkers with a 2,000-step burnin and then a 4,000-step production run for each spectral channel. We checked that the MCMC had converged by verifying that the Gelman-Rubin statistic was below 1.1 for each parameter in each spectral channel.  \nWe additionally checked the goodness-of-fit and statistical properties of our transit modeling in each spectral channel. We did so by first verifying that the average of the per-point flux uncertainties in each channel's lightcurve matched the standard deviation of the residuals to the bestfit transit model. We also com-  \nputed the Anderson-Darling statistic for each channel's lightcurve residuals to check that the residuals themselves appeared Gaussian. We did not find statistically significant non-Gaussianity in the residuals to our spectroscopic lightcurve fits.", 'B.2.2. Eureka!': "Our astrophysical model consisted of a starry (Luger et al. 2019) transit model on which we placed a broad, uninformative prior on the transit depth. We fixed our orbital parameters to those from the joint broadband fit. In addition to the linear trend in time, our systematic model included a linear correlation with the changes in the spatial position and spatial PSF-width. We also fitted a white-noise multiplier to ensure a reduced chisquared of 1 and avoid artificially constrained posteriors. No 'mirror tilt' events (Schlawin et al. 2023) were evident in either observation.  \nWe sampled the posteriors using PyMC3 's No U-Turns Sampler (Salvatier et al. 2016) with two independent chains, each taking 4,000 tuning draws and 3,000 posterior samples with a target acceptance rate of 0.85. We confirmed that the fits had converged by ensuring the Gelman-Rubin statistic was below 1.1 and by visual inspection of the two chains. Our best-fit values and uncertainties were computed using the 16th, 50th, and 84th percentiles of the PyMC3 samples. For F322W2, we find that the noise is generally ∼ 20-30% above the photon noise limit, while F444W is ∼ 2-10% above the photon noise limit. Visual examination of the Allan variance plots (Allan 1966) for both observations show no evidence for residual red noise.", 'B.2.3. tshirt': 'For our tshirt fitting, we modeled the lightcurves using a starry transit model (Luger et al. 2019). We included two linear background trends: one with respect to time and one with respect to the focal plane array housing (FPAH) temperature. For the FPAH temperature, we used the JWST engineering database to look up the IGDP NRC A T LWFPAH1 telemetry at the time of the observations. Since the telemetry data points are collected at a different cadence than our observations, we fit a 5th-order polynomial to the raw temperature values and used this to interpolate the FPAH temperatures at the time of the observations in our data. For the tshirt fitting, we binned the lightcurve and corresponding FPAH temperatures down to 300 points (roughly a factor of 5 binding for the F322W2 data and a factor of 3 binning for the F444W data) for faster lightcurve sampling.  \nSimilar to the other reductions, we fixed the orbital properties of GJ 3470 b to the values determined in the  \nFigure 7. A comparison of the three different NIRCam data reductions. The NIRCam F322W2 and F444W transmission spectra from all three of the data reduction pipelines, smoothed using a 3-point wide boxcar filter for clarity. The F322W2 spectra have all been offset down by 163 ppm (see discussion in the text). All three spectra show spectra features with similar shapes and amplitudes. Compared to our nominal Pegaus spectra, the Eureka! and tshirt results differ by χ 2 / dof = 0 . 13 for Pegasus vs. Eureka! , and χ 2 / dof = 0 . 21 for Pegasus vs. tshirt .  \n<!-- image -->  \njoint broadband fit. The tshirt fit included an additional noise-scaling parameter, which scaled the overall uncertainties on the flux measurements to allow for an increased uncertainty over the photon and read-noise expectations.  \nTo fit the spectroscopic lightcurves, we first performed a likelihood maximization to determine an initial best fit. We used this initial solution to identify and discard flux outliers that were more than 5 σ from the model. We then used the remaining points to perform a No-U Turns sampling routine via pymc3 (Salvatier et al. 2016). We sampled the lightcurves with a target acceptance of 0.90, 3,000 tuning samples and 3,000 draws from the posterior with two chains.', 'B.2.4. Comparison Between the JWST Reductions': 'The motivation for conducting three different reductions of the NIRCam data was to check the robustness of our final transmission spectrum against different choices and assumptions made during the image calibration, spectral extraction, and lightcurve fitting stages of our analyses. Generally, we find that all three reductions give the same general transmission spectrum with approximately the same depth uncertainties (Figure 2 and 7). The mean depth uncertainty across the entire NIRCam wavelength range is approximately 60 ppm at our constant wavelength ∆ λ = 0 . 015 µ m binning, which is larger than the mean point-to-point difference between the Pegasus and Eureka! spectra (28 ppm) and the Pegasus and tshirt spectra (29 ppm). A χ 2 comparison of the three spectra gives χ 2 / dof = 0 . 13 for  \nPegasus vs. Eureka! , and χ 2 / dof = 0 . 21 for Pegasus vs. tshirt . The per-point uncertainties on the F322W2 and F444W spectra from all three reductions were similar.', 'B.3. Depth Offset Between F322W2 and F444W': "There is some overlap between the NIRCam F322W2 and F444W bandpasses. In these observations, this overlap region extended from 3.89-3.95 µ m and contained four depth measurements from the F322W2 observation, and four depth measurements from the F444W observation. Visual inspection of the spectra showed a noticeable offset between the F322W2 and F444W data within this overlap region in all three of our reductions. For our nominal Pegasus results the F322W2 data in the overlap region were 163 ± 39 ppm higher than the F444W data. Both the Eureka! and tshirt spectra showed similar offsets.  \nWe used the Spitzer 3.6 µ m and 4.5 µ m transit depths we measured as a part of the joint broadband fit to determine which of the two NIRCam spectra was likely at the correct overall transit depth. To do so, we estimated a 'Spitzer equivalent' transit depth at 3.6 µ m using a portion of our F322W2 data and one at 4.5 µ m using our F444W data. We used the Spitzer/ irac filter and instrumental throughput curves provided by IPAC to account for Spitzer's varying sensitivity across each bandpass appropriately.  \nWe find that the F444W spectrum almost perfectly matches the Spitzer 4.5 µ m transit depth, while the F322W2 spectrum is significantly higher than the  \nFigure 8. Our fit to the spectral energy distribution of GJ 3470. Using catalog 2MASS and AllWISE measurements of GJ 3470's brightness along with the Gaia DR3 parallax to the system we modeled the spectral energy distribution from the star GJ 3470. We measure the stellar radius to be R ∗ = 0 . 50 ± 0 . 01 R ⊙ and its mass to be M ∗ = 0 . 44 ± 0 . 04 M ⊙ . Using the results from our lightcurve fitting, this allows us to refine the radius and mass of the planet GJ 3470 b to be R p = 4 . 22 ± 0 . 09 R ⊕ and M p = 11 . 1 ± 0 . 9 M ⊕ .  \n<!-- image -->  \nSpitzer 3.6 µ m transit depth. Specifically, the F444W spectrum gives a '4.5 µ m equivalent' depth that is lower by 17 ± 39 ppm than the measured Spitzer 4.5 µ m depth. The F322W2 spectrum gives a '3.6 µ m equivalent' depth that is higher by 144 ± 33 ppm than the measured Spitzer 3.6 µ m depth. This comparison to an independent measurement indicates that the F444W spectrum is likely at the correct overall transit depth, and the F322W2 spectrum is offset high.", 'B.4. HST Spectroscopic Transit Fitting': "To fit the spectroscopic transit depths in our rereduction of the existing WFC3 G141 data, we used a BATMAN (Kreidberg 2015) transit model and standard detrending techniques for WFC3 timeseries data (Beatty et al. 2017). As with our analysis of the JWST observations we fixed the transit center time, orbital period, scaled semi-major axis, orbital period, eccentricity, and the argument of periastron to the bestfit values from the joint broadband fit, and used a quadratic limb-darkening law with coefficients derived from the ATLAS9 (Castelli & Kurucz 2003) stellar models. We fit each spectroscopic channel individually, with the planet-to-star radius ratio and limb-darkening coefficients as the only free astrophysical parameters. We simultaneously fit the spectroscopic lightcurves from both visits for each channel, using the same value of R p /R ∗ and limb-darkening for both.  \nTo model the spacecraft and detector systematic trends present in the WFC3 data, we used an exponential ramp within each individual orbit, and a linear temporal trend across each visit, of the form  \nF detrend = ( mt V + n ) ( 1 -Ae t O τ ) . (B1)  \nHere, t O is the time since the start of each orbit's observations, t V is the time since the start of the visit, and A and τ , and m and n are fitting coefficients for the exponential ramp and linear trend respectively. In addition to fitting the systematics in this way, we also did not use the data from the first orbit within each visit nor the first exposure from each orbit when fitting the data.  \nThe measured transmission spectrum from our reanalysis of the WFC3 observations is given in Table 5.", 'C. STELLAR SED FITTING AND ABSOLUTE MASSES AND RADII': "We used a set of catalog magnitudes for the GJ 3470 system to fit an SED model to the stellar emission. In conjunction with the Gaia (Gaia Collaboration et al. 2023) parallax for the system, this allowed us to estimate a stellar radius for the star GJ 3470.  \nWe used catalog 2MASS JHK (Skrutskie et al. 2006) and AllWISE W1 and W2 (Cutri & et al. 2014) magnitudes for our SED fits. Our SED model used four different physical parameters: the stellar effective temperature, the stellar radius, the amount of visual extinction  \nto GJ 3470, and the system's parallax. We assumed log( g ) = 4 . 5 and [Fe / H] = 0 . 0 for the star GJ 3470, which are both within 0.2 dex of the surface gravity and metallicity measured by Biddle et al. (2014) In practice, the precise values of log( g ) and [Fe / H] do not significantly affect the results from the SED fit (Stevens et al. 2018).  \nWe imposed a Gaussian prior on T eff based on the spectroscopic measurements fromBiddle et al. (2014) of T eff = 3650 ± 50 K with the associated 1 σ uncertainties as the prior width. Similarly, we imposed a Gaussian prior on the parallax to the system using the Gaia DR3 (Gaia Collaboration et al. 2023) parallax of π = 34 . 0172 ± 0 . 0255 mas. For the amount of visual extinction to GJ 3470 we imposed a prior of A V = 0 . 03 ± 0 . 01, which comes from measurements (Schlafly & Finkbeiner 2011) of the excess reddening towards GJ 3470 of E ( B -V ) = 0 . 027 ± 0 . 01, and assuming R V = 3 . 1.  \nTo model the SED, we used BHAC15 spectra (Baraffe et al. 2015) for the stellar SED. We computed a grid of surface luminosity magnitudes, corresponding to the bandpasses of the catalog magnitudes, for a range of T eff values. Since the BHAC15 models step by 200K in T eff we used cubic spline interpolation to estimate model magnitudes in between the points provided by the model atmospheres. We then scaled the interpolated surface magnitudes for the star by R ∗ /d - where d is the distance to the star - to determine the apparent bolometric flux of the SED at Earth. We then applied a simple R = 3 . 1 extinction law scaled from the value of A V , to determine the extincted bolometric flux of the SED model.  \nThis stellar SED fitting (Figure 8) allowed us to measure the radius of the star GJ 3470 to be R ∗ = 0 . 50 ± 0 . 01 R ⊙ . We then used the mean value of R P /R ∗ as measured in our joint broadband transit fits to estimate the planetary radius of GJ 3470 b to be R p = 4 . 22 ± 0 . 09 R ⊕ . Using the stellar density we measure as a part of the joint broadband transit fits, we constrain GJ 3470's stellar mass to be M ∗ = 0 . 44 ± 0 . 04 M ⊙ , and using a mass-ratio measured via radial velocity observations (Kosiarek et al. 2019), we estimate a planetary mass of M p = 11 . 1 ± 0 . 9 M ⊕ . Our atmospheric retrievals use these absolute stellar and planetary radii and masses.", 'D. STELLAR VARIABILITY': "Previous photometric monitoring of GJ 3470 during the 2012-2013 season showed that the star itself is slightly variable at the presumed stellar rotation period of 21 days (Biddle et al. 2014). This same monitoring campaign now spans from 2012 to 2023 and has given us  \na more precise measurement of stellar variability (Figure 9).  \nWe observed the GJ 3470 system with the Tennessee State University Celestron 14-inch (C14) automated imaging telescope (AIT) located at Fairborn Observatory in southern Arizona (Eaton et al. 2003). From 2012 to 2023 we collected 1376 good observations using a Cousins R filter with an SBIG STL-1001E CCD camera. Each nightly observation consisted of three to ten consecutive exposures of the GJ 3470 field of view, which we use to perform differential aperture photometry on GJ 3470. Further details of the observing and data reduction procedures have been described previously in Sing et al. (2015).  \nUsing these new additional observations, we find that the GJ 3470 system shows sinusoidal brightness modulations with an amplitude of 2545 ± 130 ppm at a period of P = 21 . 6239 ± 0 . 0014 days. This variability amplitude is roughly half the value measured previously using only the 2012-2013 data (Biddle et al. 2014) and roughly half the variability amplitude used in previous modeling of how stellar variability could affect GJ 3470 b's transmission spectrum in Benneke et al. (2019b). Note that we have normalized the seasonal means of the photometric observations to all equal unity. This allows us to estimate the short-term (tens of days) variability, but not multi-year variability, of GJ 3470.  \nIf we conservatively assume that the amplitude of GJ 3470's rotational variability is constant as a function of wavelength, then this will introduce at most a 30 ppm offset between our various observations. In reality, the amplitude of stellar variability decreases significantly at longer wavelengths. Benneke et al. (2019b) modeled the wavelength dependence of GJ 3470's variability and its effect on transmission spectra using a variability amplitude approximately twice that which we measured, and found that it was at most 20 ppm in the WFC3 bandpass and < 10 ppm over the F322W2 and F444W wavelength range. Since even our conservative estimate for the effect of stellar variability - assuming a wavelengthindependent amplitude - is on the order of (or smaller) than our estimated depth uncertainties and the interinstrument offsets, we do not expect the variability seen in the GJ 3470 system to affect our results significantly.", 'E. ATMOSPHERIC MODELING AND RETRIEVALS': 'In order to derive quantitative atmospheric constraints we perform a series of atmospheric retrievalsBayesian parameter estimation with a planetary atmosphere/spectral model. Our setup/philosophy follows similar to that described in Bell et al. (2023)  \nFigure 9. Photometric monitoring of the GJ 3470 system from 2012 to 2023 shows sinusoidal brightness modulations with an amplitude of 2545 ± 130 ppm at a period of P = 21 . 6239 ± 0 . 0014 days in a Cousins R filter. The amplitude of stellar variability decreases significantly at longer wavelengthsBenneke et al. (2019b) and we expect this variability to affect our measured transit depths at < 10 ppm over the F322W2 and F444W wavelength range.  \n<!-- image -->  \nand Welbanks et al. submitted. Briefly, we perform both 1D-radiative-convective-photochemical equilibrium grid-based (1D-RCPE) retrievals and more flexible, classic, free-retrievals, each summarized below (Welbanks & Madhusudhan 2021).', 'E.1. 1D-RCPE Grid-Retrieval': 'We use the ScCHIMERA (Piskorz et al. 2018; Mansfield et al. 2021; Bell et al. 2023) 1D-radiative-convectivephotochemical-equilibrium solver (1D-RCPE) to produce a physically self-consistent grid of atmospheric models that are then fit to the spectra using the methods outlined in Bell et al. (2023). The grid is a function of metallicity (0 ≤ [M/H] ≤ 3.0 in steps of 0.125), C/O (0.01 ≤ C/O ≤ 0.6, using 8 spacings), internal temperature (100 ≤ T int ≤ 350 in steps of 50 K), and a vertically constant eddy mixing (7 ≤ log 10 K zz ≤ 10 in steps of 0.5), at a fixed effective irradiation temperature (at the planetary equilibrium temperature-653 K). These resulted in 8400 converged 1D-RCPE model atmospheres. For the photochemistry/kinetics (using the VULCAN (Tsai et al. 2021) chemical kinetics code), we use the H-O-CN-S reaction list (Tsai et al. 2022) and the MUSCLES (France et al. 2016; Youngblood et al. 2016; Loyd et al. 2016) GJ 176 spectrum as a proxy for GJ 3470.  \nAs in previous work (Bell et al. 2023), we fit for the metallicity, C/O, internal temperature, eddy diffusivity, planetary radius and planetary radius reference pressure, gray and power law patchy clouds (Welbanks & Madhusudhan 2021), and stellar heterogeneity (Iyer & Line 2020). Parameter estimation is per-  \nformed with the pymultinest (Buchner et al. 2014) Nested Sampling tool. Within the Nested Sampling, the pre-computed 1D-RCPE model atmospheres (temperature, molecular weight, and gas volume mixing ratio profiles) are interpolated (using the python/scipy RegularGridInterpolator) to the sampled metallicity, C/O, internal temperature, eddy diffusivity. The interpolated model atmosphere at each step/sample-along with the current cloud, reference radius/pressure, and stellar heterogeneity parameters are used to compute a 1D highresolution (R=100K) transmission spectrum. This spectrum is then binned down to the data resolution to compute the chi-square log-likelihood used by the sampler. In summary, these grid-based constraints suggest a metallicity 125 ± 40 times Solar, a slightly sub-solar C / O = 0 . 35 ± 0 . 1, weak vertical mixing, and a poor constraint on the internal temperature. A power law scattering slope is needed to explain the join the WFC3to-NIRCam continuum/slope. Stellar heterogeneity also appears to play little role in the shape of the spectrum as the spot-covering fraction pushes towards the lower prior bound of zero. Main text Figure 4 shows the agreement between the molecular abundances retrieved from the free retrieval and those predicted by the retrieved grid-based retrieval parameters suggesting that the inferred compositional state is plausible.', 'E.2. Free-Retrievals': "We also explored GJ 3470 b's atmospheric properties using more flexible and agnostic forward models in a Bayesian inference procedure. These models, known  \nas free-retrievals, are methods to retrieve the chemical abundances of different gases, the vertical temperature structure, and the cloud/haze properties on the planetary atmosphere. Compared to the 1D-RCPE GridRetrieval explained above, these methods do not assume any physico-chemical equilibrium conditions, and instead aim to capture the atmospheric conditions directly through a series of parameters for the chemistry and physical conditions of the atmosphere without expectations of physical consistency. These more flexible approaches provide the opportunity to capture conditions that otherwise would be prohibited by selfconsistent models, such as combinations of gases not considered under chemical equilibrium. Nonetheless, caution must be exercised in the interpretation of freeretrievals as model assumptions may contribute to biased and unphysical atmospheric estimates (Welbanks & Madhusudhan 2022).  \nWe employ the retrieval framework Aurora (Welbanks & Madhusudhan 2021) capable of modeling and interpreting transmission and emission spectra of transiting exoplanets (Welbanks & Madhusudhan 2022; Bell et al. 2023). A detailed description of Aurora's transmission modeling approach is described in Welbanks & Madhusudhan (2021). Briefly, the atmospheric model solves radiative transfer for a parallel-plane atmosphere under hydrostatic equilibrium for a transmission geometry. The atmospheric model considers a one-dimensional model atmosphere spanning from 10 -9 bar to 100 bar, divided into 100 layers uniformly spaced in logarithmic pressure space. The vertical temperature structure is parameterized following Madhusudhan & Seager (2009). The model simultaneously retrieves the reference pressure and reference radius for the planetary gravity of log( g p ) = 2 . 78 cgs (Table 6).  \nThe atmospheric models assume uniform mixing ratios for H 2 O, CH 4 , NH 3 , CO, CO 2 , SO 2 , and HCN using independent free parameters for each gas' volume mixing ratio. The presence of clouds and hazes is incorporated via a one-sector parameterization (Welbanks & Madhusudhan 2021), which treats the combined spectroscopic effect of an optically thick cloud deck at a pressure P cloud together with hazes following an enhancement to Rayleigh-scattering (Lecavelier Des Etangs et al. 2008) in a linear combination with a cloud-free atmosphere (Line & Parmentier 2016).  \nThe Bayesian inference is performed using nested sampling (Skilling 2004; Feroz et al. 2009) via PyMultiNest (Buchner et al. 2014) using 500 live points in the sampling with an evidence tolerance of 0.5 and a sampling efficiency of 0.8 . Each forward model in the sampling was calculated using line-by-line opacity sam-  \npling at a spectral resolution of 20,000 and then convolved and binned to the resolution of the observations. In total, the sampling is performed over 19 parameters: 7 molecular gases, 6 for the pressure-temperature structure of the planet, 1 for the reference pressure, 1 for the reference radius, and 4 for the presence of inhomogeneous clouds and hazes.  \nFigure 3 shows the retrieved transmission spectrum and the contributions of the detected gases in GJ 3470 b's atmosphere. The retrieved volume mixing ratios from the free-retrieval are consistent with the gas profiles from the 1D-RCPE inferences, as shown in Figure 4.  \nBoth free retrievals confirm the detections of H 2 O, CH 4 , CO 2 , and SO 2 with detection significances of ∼ 4 σ or greater. We note that the term 'detection significance' refer to a model preference between a reference model considering all gases and nested models for which each individual species is removed in turn (Benneke & Seager 2013; Welbanks & Madhusudhan 2021). As such, any quoted detection significance is dependent on the choice of models and choice of model priors. For GJ 3470 b, the agreement in retrieved molecular abundances and confirmation of the strong evidence for the detection of these species regardless of model assumptions (e.g., free abundances vs. radiative convective-photochemicalequilibrium) confirms the robustness of the interpretation of the planet's spectrum. Table 1 shows the retrieved parameters, priors, and detection significances where applicable.  \nTable 2. Median Values and 68% Confidence Intervals for the Joint Broadband Fit  \nb The broadband WFC3, F322W2, and F444W data were included in the joint broadband fitting process, but we used the transmission spectra from these instruments in our atmospheric modeling.  \nTable 3. GJ 3470 b's NIRCam/F322W2 Transmission Spectrum  \nThe full machine-readable version of this table is available online.  \nNote: These values have been offset by -163 ppm as discussed in Section 2.4 and Appendix B.3.  \nTable 4. GJ 3470 b's NIRCam/F444W Transmission Spectrum  \nThe full machine-readable version of this table is available online.  \nTable 5. GJ 3470 b's HST/WFC3 Transmission Spectrum  \nTable 6. Median Values and 68% Confidence Intervals for the Stellar and Planetary Properties"}