JWST/NIRCam Transmission Spectroscopy of the Nearby Sub-Earth GJ 341b

We present a JWST/Near Infrared Camera (NIRCam) transmission spectrum from 3.9 to 5.0 μm of the recently validated sub-Earth GJ 341b (R P = 0.92 R ⊕, T eq = 540 K) orbiting a nearby bright M1 star (d = 10.4 pc, K mag = 5.6). We use three independent pipelines to reduce the data from the three JWST visits and perform several tests to check for the significance of an atmosphere. Overall, our analysis does not uncover evidence of an atmosphere. Our null hypothesis tests find that none of our pipelines’ transmission spectra can rule out a flat line, although there is weak evidence for a Gaussian feature in two spectra from different pipelines (at 2.3 and 2.9σ). However, the candidate features are seen at different wavelengths (4.3 μm versus 4.7 μm), and our retrieval analysis finds that different gas species can explain these features in the two reductions (CO2 at 3.1σ compared to O3 at 2.9σ), suggesting that they are not real astrophysical signals. Our forward-model analysis rules out a low-mean-molecular-weight atmosphere (<350× solar metallicity) to at least 3σ, and disfavors CH4-dominated atmospheres at 1–3σ, depending on the reduction. Instead, the forward models find our transmission spectra are consistent with no atmosphere, a hazy atmosphere, or an atmosphere containing a species that does not have prominent molecular bands across the NIRCam/F444W bandpass, such as a water-dominated atmosphere. Our results demonstrate the unequivocal need for two or more transit observations analyzed with multiple reduction pipelines, alongside rigorous statistical tests, to determine the robustness of molecular detections for small exoplanet atmospheres.

The M-dwarf opportunity is one of the most promising routes to measure the atmospheres of terrestrial, and potentially habitable, exoplanets (e.g.Charbonneau & Deming 2007).This opportunity arises from the favorable planet-to-star radius ratios, close-in habitable zones, and frequent transits of rocky planets transiting M dwarfs.However, the M-dwarf opportunity must contend with the high and prolonged X-ray and EUV emission from M-dwarf host stars (e.g., Pizzolato et al. 2003;Preibisch & Feigelson 2005;Shkolnik & Barman 2014;Peacock et al. 2019), which could cause significant exoplanet atmospheric loss (e.g., Lammer et al. 2003;Owen & Jackson 2012;Owen & Mohanty 2016).
In order to assess the M-dwarf opportunity, it is necessary to determine which, if any, terrestrial planets around M dwarfs are capable of retaining an atmosphere.One way to address this question is to determine whether M dwarf exoplanets are separated by a 'cosmic shoreline', a division in insolation-escape velocity space (I-v esc ) separating bodies with atmospheres from those without (Zahnle & Catling 2017).In the Solar System, Zahnle & Catling (2017) demonstrated that this dividing line follows an I ∝ v esc 4 relation.Following the successful launch and commissioning of JWST (Gardner et al. 2023;McElwain et al. 2023;Rigby et al. 2023), we are now able to expand the I-v esc parameter space over which we can address the existence of a cosmic shoreline to terrestrial planets around M-dwarf stars.This is the focus of JWST program 1981 (PIs: K. Stevenson and J. Lustig-Yaeger), which is obtaining JWST transmission spectra of five exoplanets around M dwarfs to determine whether a cosmic shoreline exists for M-dwarf systems.
Our first study (Lustig-Yaeger & Fu et al. 2023b) ruled out Earth-like, hydrogen/helium, water-and methane-dominated clear atmospheres for the Earthsized exoplanet LHS 475b.This is consistent with predictions from the cosmic shoreline.This study also confirmed LHS 475b as a bona fide exoplanet, promoting it from a planet candidate in parallel with the efforts of Ment et al. (2023).
Our second study (Moran & Stevenson et al. 2023) was of the super-Earth GJ 486b which sits close to the shoreline predicted by Zahnle & Catling (2017) and thus provides an important test of its existence.We tentatively detected water in GJ 486b's transmission spectrum, however, our NIRSpec/G395H observations could not determine the origin of the water, be it in the planet's atmosphere or the stellar atmosphere.
Most recently, in May & MacDonald et al. 2023, we observed two transits of the super-Earth GJ 1132b, another planet near the shoreline.Our study revealed the importance of multi-visit observations for small transiting exoplanets, since one visit was suggestive of spectral features while the other produced a flat transmission spectrum.Given the inconsistent transmission spectra between our two visits, we were unable to make a conclusive determination of the planet's likely atmospheric composition, although neither visit was consistent with a previous report of HCN (Swain et al. 2021).
Here, we extend our program to GJ 341b (TOI 741.01), a planet that was recently validated by DiTomasso et al. (2023) and which orbits a nearby bright M1 star (K = 5.6, 10.4 pc, Gaia Collaboration et al. 2016, 2023).The TESS-derived planetary radius is 0.92±0.05R ⊕ , making this a sub-Earth, with a 3σ upper limit on the planet's mass of 4.5 M ⊕ (DiTomasso et al. 2023).The mass-radius relation of Chen & Kipping (2017) gives a planet mass of 0.72 ± 0.14 M ⊕ .Assuming a Bond albedo of zero and uniform heat redistribution, the planet's equilibrium temperature is 540 K, while its irradiation temperature is 760 K. Given the planet's escape velocity of 10 km s −1 , the cosmic shoreline relation predicts that this planet has not retained a substantial atmosphere.Our study seeks to test this prediction.
The paper is laid out as follows: in Section 2 we describe the observational setup.In Section 3 we explain our data reduction and light curve fitting, using three independent pipelines.We present our transmission spectra, flat line rejection tests, atmospheric forward models and atmospheric retrievals in Section 4. We discuss our results and present our conclusions in Section 5.

OBSERVATIONS
We observed three transits of GJ 341b across three visits with JWST's Near Infrared Camera (NIRCam, Rieke et al. 2023) as part of GO program 1981 (PIs: K. Stevenson and J. Lustig-Yaeger).Unlike the other GO 1981 targets that were observed with NIRSpec, NIRCam was selected for GJ 341 because this bright star (K mag = 5.6) saturates NIRSpec in one group.The three visits occurred on March 2, March 10, and April 17, 2023 UT, with each visit lasting 5.25 hours.We used the F444W filter, covering a wavelength range from 3.8-5.1 µm at a native spectral resolution of R ≈ 1650, and the SUB-GRISM128 subarray, giving us 2048 pixel columns and 128 pixel rows.Given the brightness of our target, we used 3 groups per integration and acquired 4652 integrations during each of our three visits.The effective integration time was 3.38 s.

DATA REDUCTION
As with the previous studies in our program (Lustig-Yaeger & Fu et al. 2023b;Moran & Stevenson et al. 2023;May & MacDonald et al. 2023), we used three different pipelines to reduce our data independently: Eureka!, Tiberius, and Tswift.We describe the approach taken with each pipeline in more detail below and include a table in the appendix that summarizes each pipeline's approach (Table 6).

Eureka!
The Eureka! data reduction pipeline (Bell et al. 2022) is well established, with published JWST analyses demonstrating reliable results from numerous observations (The JWST Transiting Exoplanet Community Early Release Science Team et al. 2022;Ahrer et al. 2023;Alderson et al. 2023;Rustamkulov et al. 2023;Lustig-Yaeger & Fu et al. 2023b;Moran & Stevenson et al. 2023).As described below, we adopt similar methods for the analysis of GJ 341b.
Starting from the uncal.fitsfiles, we masked a region within 30 pixels of the measured trace from pixel column 640 to 2048.We then performed a doubleiteration, group-level background subtraction by first fitting and removing the mean flux on a column-by-column basis and, second, fitting and removing the mean flux on a row-by-row basis.Unlike NIRSpec and NIRISS, NIRCam detectors exhibit 1/f noise along pixel rows, thus necessitating this additional row-by-row step.This orientation of the 1/f noise is because NIRCam's rapid readout direction aligns with the wavelength dispersion direction (Schlawin et al. 2020).In Stage 1, we skipped the jump detection step because our observations only had 3 groups per integration, which makes the jump step unreliable.We processed the rateints.fitsfiles through the regular jwst Stage 2 pipeline.
In Stage 3, Eureka!performed a full-frame outlier search with a 4σ rejection threshold.We then computed a clean median frame and corrected the spectrum's curvature by rolling pixel columns to bring the trace to the center of the detector.We skipped any additional background subtraction and applied optimal spectral extraction (Horne 1986) using an aperture that was within 3 pixels of the measured trace.We tested a range of aperture sizes and found that a 3-pixel half-width minimized the white light curve residuals.To calculate our wavelength solution, we used a polynomial from Flight Program 1076 (E.Schlawin, private communication).
We extracted the flux from 3.9-5.0µm, splitting the light into 55 spectroscopic channels, each 0.02 µm in width.For each light curve, we applied a box-car filter of length 50 integrations and flagged any 3.5σ outliers.This procedure removed < 2% of all integrations for a given light curve.
We fitted all three white light curves simultaneously.We used batman (Kreidberg 2015) to fit the system parameters and fixed the quadratic limb-darkening coefficients to those provided by ExoTiC-LD (Grant & Wakeford 2022) from 3D stellar models (Magic et al. 2015).For our stellar parameters, we used T eff = 3745 K and log g = 4.71 g cm −2 and assumed a solar metallicity ([Fe/H] = 0) (ExoMAST).Across the three visits, we fitted for a single planet-star radius ratio (R P /R * ), orbital inclination (i), semi-major axis (a/R * ), and transit time (T 0 ).The period (P) was held fixed to 7.5768707 days (ExoMAST).For each visit, we fitted an independent exponential times linear ramp model to account for light curve systematics.We also renormalized the error bars from each visit in order to achieve a reduced chisquared of unity.In total, there were 19 free parameters.Table 1 lists our best-fit system parameters.
When fitting the spectroscopic light curves, we fixed the orbital inclination, semi-major axis, and transit time to the white light curve best-fit values.Again, we fixed the quadratic limb-darkening parameters to the values provided by ExoTiC-LD for each spectroscopic channel.The planet-to-star radius ratio was the only physical free parameter.For each spectroscopic channel, we fitted an independent exponential times linear ramp model to account for light curve systematics.Including the error bar renormalization term, each spectroscopic fit had six free parameters.For both the white and spectroscopic light curves, we clipped the first 300 integrations from each visit to avoid fitting the exponential decrease in flux.Figure 1 shows our fits to the Eureka!white and spectroscopic light curves, with the Allan variance plot of the residuals shown in Figure 2.

Tiberius
The Tiberius reduction began with the uncalibrated fits files (uncal.fits).We ran these through the stage 1 detector processing stages of the jwst pipeline.However, we did not perform the jump step as this has been found to increase the noise in time series spectrophotometry (e.g., Rustamkulov et al. 2023).Furthermore, unlike in the application of Tiberius to NIRSpec data (e.g., Rustamkulov et al. 2023;Lustig-Yaeger & Fu et al. 2023b;Moran & Stevenson et al. 2023) and in our Eureka!reductions of the data presented here, we did not apply a 1/f correction to the data.This is because  1. GJ 341b's system parameters from the three independent reductions' white light curve fits.We also include the parameters from the independent TESS and radial velocity analysis (DiTomasso et al. 2023).
Figure 1.The light curves and best fit models resulting from the Eureka!pipeline.The columns show, from left to right: transit 1, transit 2, transit 3. The rows show, from top to bottom: 1) the white light curve and residuals, 2) the spectroscopic light curves, 3) the best-fit models to the spectroscopic light curves, 4) the residuals from the spectroscopic light curve fitting.
the 1/f noise is parallel to the dispersion dimension for NIRCam data, complicating its removal.Following the stage 1 detector processing, we worked with the integration level images (gainscalestep.fitsfiles).
While these were processed through the gain scale step, the jwst pipeline does not apply a gain correction if the default gain setting is used.For this reason, we worked in units of DN/s throughout the rest of our analysis.Since we work with normalized light curves, and rescale our photometric uncertainties, the gain correction is not necessary.
Before extracting the stellar flux from each image, each pixel was compared to the median value for that pixel across the time series and > 5σ outliers were replaced by the median for that pixel.We then oversampled each pixel by a factor of 10 in the spatial dimension using a flux-conserving interpolation, which allows us to more precisely place the aperture boundaries.We found this step to reduce white noise by 10 % in previous tests on JWST data (Rustamkulov et al. 2023).Next, we traced the stellar spectrum using a fourth order polynomial between columns 700 and 2042 (inclusive, zeroindexed) and performed standard aperture extraction.We subtracted the background at every column using the median value across the full column, after masking a total of 45 pixels centered on the stellar trace, and summed up the stellar flux in an aperture with a full width of 5 pixels.We experimented with wider apertures and found these led to consistent transmission spectra and light curves but with greater noise.We fitted our white light curves from the three visits simultaneously after subtracting one orbital period (7.576863 days, DiTomasso et al. 2023) from transit 2's time array and six orbital periods from transit 3's time array so that all time arrays could be fitted with a single T 0 , which corresponded to transit 1.The 16 free parameters in our white light curve fits were T 0 , i, a/R * , R P /R * and the 12 parameters defining our systematics model.We defined this model as a single component exponential ramp multiplied by a linear-in-time polynomial for each visit.The period was held fixed to 7.576863 days (DiTomasso et al. 2023).We used batman (Kreidberg 2015) to generate our analytic transit light curves.
For our limb darkening calculations, we adopted T eff = 3770 ± 40 K and log g = 4.72 ± 0.02 g cm −2 (DiTomasso et al. 2023) and assumed a solar metallicity ([Fe/H] = 0).We used a quadratic limb darkening law with coefficients held fixed to values calculated from 3D stellar atmosphere models (Magic et al. 2015) using ExoTiC-LD (Grant & Wakeford 2022) and our own custom instrument throughput.This was constructed by dividing our observed stellar spectrum by the appropriate Phoenix model spectrum (Husser et al. 2013) and smoothing the result with a running median, following the approach of Alderson et al. (2023).
We explored the parameter space with a Markov chain Monte Carlo (MCMC) using the emcee Python package (Foreman-Mackey et al. 2013).We ran three sets of chains each with 50,000 steps and 200 walkers.After the first run, we discarded the first 25,000 steps and calculated the median values from the second 25,000 steps of our chains which were used as starting points for the second run.After the second run, we again discarded the first 25,000 steps and used the best-fitting model to rescale our photometric uncertainties to give χ 2 ν = 1.Our third and final run used the rescaled photometric uncertainties.The medians, 16th and 84th percentiles of our posteriors from this third run are reported in Table 1.These posteriors were constructed after discarding the first 100 steps and thinning by a factor of 10.
Our spectroscopic light curves were created using 57 20-pixel-wide wavelength bins between 3.906 and 5.004 µm and we adopted the same wavelength solution as used in the Eureka!reduction.For our spectroscopic light curve fits, we fixed the system parameters (T 0 , i and a/R * ) to the Tiberius values in Table 1 and fitted each visit independently, not jointly.Prior to fitting, we divided each spectroscopic light curve by a visit-dependent common noise model.This noise model was defined as WLC v /tm v where WLC v is the visit-dependent white light curve and tm v is the visitdependent best-fitting transit (batman) model.For our spectroscopic light curves, the exponential coefficients did not converge for every wavelength bin across each visit, despite experimenting with different priors.For this reason, at the spectroscopic light curve stage, we clipped the first 600 integrations (40 minutes) from each visit, after removing the common noise, and fitted the light curves with a linear polynomial and no exponential.This gave comparable residual red noise to the fits where the exponential coefficients did converge, indicating that the ramp was adequately removed by clipping the first 600 integrations.Figure 11 in the Appendix shows the fits to our white and spectroscopic light curves, with the Allan variance plot of the residual red noise from our white light curve fits shown in Figure 2.

Tswift
Transit swift (Tswift) is a JWST data reduction routine developed to swiftly analyze time-series observation of transiting exoplanet atmospheres.Starting from the uncal.fitsimages, we used the default jwst pipeline with steps.jumpturned off to produce stage one rateints.fitsimages.Next, we performed rowby-row and column-by-column background subtraction (Schlawin et al. 2020) for each integration by subtracting out the median of the non-illuminated regions from each row.Next, we cross-correlated the spectral line spread function with each column to determine the vertical position of the spectral trace for that column.We then used a third-order polynomial to fit the vertical position versus the column number, determining the spectral trace location for each integration.The spectrum was extracted using a full-width 6-pixel wide aperture centered on this polynomial fit as it minimizes the scatter in the white light curve.Given the brief exposure time for each integration and the high number of integrations per visit, cosmic rays or hot pixels have a minimal impact on the transit light curve.Furthermore, we identified and discarded any outliers in time exceeding 5σ based on the transit light curve baseline scatter for each column.Finally, we combined all wavelength channels to produce the white light curve of each transit.
The fitting of each individual white light curve employs a combination of batman (Kreidberg 2015) and emcee (Foreman-Mackey et al. 2013).This approach incorporates eight parameters: T 0 , R P /R * , a/R * , i, a linear time slope, a constant scaling factor, and two coefficients for the exponential ramp.The averaged best-fit parameters from the three white light curves are listed in Table 1.To determine the quadratic limb darkening parameters, we used the 3D stellar models (Magic et al. 2015) from the ExoTiC-LD (Grant & Wakeford 2022) package and subsequently fixed both quadratic terms, using the same stellar parameters as in the Eureka! and Tiberius analyses.We tested fitting for the second quadratic term while fixing the first one, and it had no effect on the final transmission spectrum.
The NIRCam detector has four amplifiers and the spectrum is dispersed across the second to fourth amplifier in the F444W observation mode (Schlawin et al. 2020).As already discussed, there is 1/f noise along the dispersion direction in our NIRCam data.In an ideal scenario, this noise could be subtracted if every row in each amplifier had unilluminated pixels.However, due to the continuous dispersion of spectra across various amplifiers, it is not feasible to correct the entire spectrum.As a result, we observed lingering 1/f noise residuals that showed correlation across wavelengths in the light curve fits.We identified that these residual patterns are predominantly consistent across all amplifiers, with a few minor patterns unique to individual amplifiers.To better remove this correlated noise, we performed a white light fit within each amplifier, and the residuals were then used as a common-mode correction during the spectroscopic light curve fits.
For our spectroscopic light curve fits, we used batman and scipy.optimize.curvefit.The model consists of five parameters: R P /R * , a constant scaling factor, two coefficients for the exponential ramp, and a scaling factor for the common-mode correction.We held T 0 , a/R * , and i fixed to the corresponding averaged best-fit values from the white light curves (Table 1).Figure 12 shows our fits to the Tswift white and spectroscopic light curves, with the Allan variance plot of the residuals shown in Figure 2.

Engineering parameters
As can be seen in the Allan variance plots from our reductions (Figure 2) there is residual red noise which becomes apparent at time scales longer than 1 minute.In order to remove this red noise, we explored detrending our light curves with engineering telemetry data timeseries from the JWST Engineering Database.The telemetry data are identified by mnemonics, which follow the naming convention X Y Z, where X corresponds to the system (I for Integrated Science Instrument Module, S for spacecraft), Y corresponds to the subsystem, instrument, controller, or component (1-4 letters), and Z corresponds to the telemetry name (2-20 characters).Throughout the rest of this section, we refer to these timeseries telemetry data sets as 'mnemonics'.
We began by whittling down the 6176 NIRCamrelated mnemonics to the 596 mnemonics of timevarying floats (as opposed to booleans and non-varying constants).We then resampled each mnemonic onto the same time-sampling as each of our three visits and proceeded to sequentially incorporate each of the 596 mnemnonics into our white light curve systematics model, for the Tiberius reduction, to test whether the addition of a mnemnonics model decreased the residual red noise.
For these fits we used a Levenberg-Marquadt alogorithm, with the systematics model being the summation of a linear polynomial describing the mnemonic and a linear-in-time polynomial, after clipping the first 600 integrations.We then compared the residual red noise in each case to our fiducial model: a linear-in-time polynomial.We defined the residual red noise as being the RMS of the residuals in bins of 200 integrations divided by the white noise expectation (equivalent to dividing the color lines in our Allan variance plots by the dashed black lines at x = 14.5 mins, Figure 2).Of the 596 possible mnemonics, we found that 112 mnemonics led to a decrease in the residual red noise compared to the fidu-cial model for all three visits.The ten mnemonics that led to the largest decrease in the RMS are given in Table 2.As shown in Table 2, these mnemonics mostly correspond to temperature.
Following this test, we ran a second test where we increased the number of mnemonics used to fit the white light curves, going from one to ten mnemonics in the order presented in Table 2.As can be seen in Figure 3, including additional mnemonics does not necessarily lead to an improvement in the residual red noise.This is likely due to some mnemonics representing similar telemetry data (e.g., INRC FA PSC LVCBRD TEMP and INRC FA PSC SCRBRD TEMP).This figure also demonstrates that the use of a mnemonic significantly improves the residual red noise in visit 2 while the changes in visits 1 and 3 are more subtle.Additionally, a combination of 10 mnemonics does not eradicate the residual red noise.As a further test, we fitted all 596 mnemonics jointly with a linear-in-time polynomial.In this case, the residual red noise was below the white noise expectations.We interpreted this as evidence that including a large number of mnemonics can lead to overfitting.
While independently fitting each visit's light curve with mnemonics demonstrated the residual red noise could be improved upon, we ended up not using the mnemonic data in our final white light curve fits.This was because we could not obtain acceptable fits when we tried to fit the three visits simultaneously with mnemonics, likely due to the fact that the amplitude of the systematics in our data is comparable to the transit depth.We advocate for a further exploration of mnemonics for exoplanets with deeper transits where the transit and systematics models are less degenerate.

Short wavelength photometry
Since our observations were conducted with NIRCam, we additionally have access to the short wavelength (SW) channel between 0.6-2.3µm.This channel acquired photometry simultaneously with the long wavelength spectroscopy described above.We extracted the short wavelength photometry using Tiberius to determine whether we could use these additional white light curves to improve the precision of our system parameters.
As with the long wavelength data, we began by processing the uncal.fitsfiles through stage 1 of the jwst pipeline.We then extracted the time-series photometry using Tiberius with an aperture 100 pixels wide and a center fixed at a pixel column of 1795.We also subtracted the background at every pixel row, which was calculated as the median of two 50-pixel-wide regions separated by 100 pixels from each side of the target aperture.The resulting SW light curves for each visit are shown in Figure 4. Given the significant noise in these light curves, we did not use these in our parameter estimation.

Measuring flux-calibrated spectra
Following our study of GJ 486b (Moran & Stevenson et al. 2023), we extracted the flux calibrated stellar spectrum of GJ 341.We did this in an attempt to constrain the covering fractions of star spots and faculae and hence determine the impact of stellar contamination on our transmission spectra.To do this, we processed our stellar spectra through stage 2 of the jwst pipeline.Specifically, taking the gainscalestep.fitsfiles, we processed these through the assign wcs, flat field, extract 2d, srctype, pathloss and photom steps.We did not perform the background step, and instead used an alternate version of our stage 1 reduction where we did perform a group-level background subtraction.
We then used Tiberius to perform standard aperture photometry on our photomstep.fitsfiles using an aperture with a full width of 20 pixels.Despite running the srctype command with the --source type POINT argument defined, the resulting units of our spectra were MJy/sr.We therefore had to multiply our spectra by the solid angle of a pixel, listed as PIXAR SR in the photomstep.fitsheader, to convert to MJy.We defined the flux uncertainties in the flux calibrated spectra as being the standard deviation of each pixel's flux across the time series.We found no significant variability in the flux-calibrated spectra from visit-to-visit (≤ 0.1 %) and within each visit (0.24 %).
We then computed multi-component forward models using the Allard et al. (2012) PHOENIX stellar models to quantify star spot and faculae covering fractions of the stellar surface for each of the three visits.We employed a weighted linear combination of three PHOENIX models to represent the background photosphere, spots (T eff ≤ T eff,photosphere -100 K), and faculae (T eff ≥ T eff,photosphere + 100 K).We assumed that all spots have the same T eff , log(g), and metallicity, as do the facu-lae, and each feature is required to not exceed 45% of the total surface.The grid of PHOENIX models used for our analysis covered T eff = 2800 -4000 K, log(g) = 4 -5.5 cm s −2 , [Fe/H] = -0.5 -0.5 providing extensive coverage of possible spot and faculae temperatures for GJ 341 assuming photospheric values near the literature quoted T eff = 3770 ± 40 K and log(g) = 4.72 ± 0.02 cm s −2 (DiTomasso et al. 2023).The models were scaled by R 2 * /dist 2 using literature values for GJ 341: R * =0.506 R ⊙ (DiTomasso et al. 2023) and d = 10.45 pc (Gaia Collaboration et al. 2021).We also smoothed and interpolated the models to be the same resolution as the observations before calculating a reduced-χ 2 .In our reduced-χ 2 calculations, we considered 868 wavelength points for each visit and eight fitted parameters (the T eff , log(g), and [Fe/H] of the photosphere, the T eff and coverage fraction of both spots and faculae, and a scaling factor).The scaling factor was multiplied by the R 2 * /d 2 term to account for uncertainty in either measured quantity and varied from 0.9 to 1.3.
Unlike in our application of this procedure to the stellar spectra of GJ 486 (Lustig-Yaeger & Fu et al. 2023b) and GJ 1132 (May & MacDonald et al. 2023), which used NIRSpec/G395H, we could not obtain an acceptable fit to the flux-calibrated NIRCam/F444W spectra (χ 2 ν ∼ 60), as shown in Figure 5.We do not know the reason for these poor fits, however, the shape of our extracted spectrum changes significantly when we apply the photometric calibration (photom) step in the jwst pipeline, suggesting that the NIRCam PHOTOM reference file may be responsible.Given the poor fits in Figure 5, and lack of contamination evident in our transmission spectra (see section 4), we do interpret the results from our analysis of the flux calibrated spectrum as evidence for stellar contamination in our data.

GJ 341b's TRANSMISSION SPECTRUM
The transmission spectra we measure for GJ 341b from the three different reductions are shown in Figure 6.This figure demonstrates that the three independent reductions produce consistent transmission spectra.For this figure, we offset the Eureka!spectrum by -15 ppm so that its median transit depth matches that of Tiberius and Tswift between 4.0 and 4.8 µm.We include the transmission spectra measured for each visit by each reduction in the Appendix (Figure 13) which demonstrates that the spectra are also consistent between visits for each reduction.In the following subsections, we further investigate the transmission spectrum of GJ 341b with flat line rejection tests, atmosphere forward models, and retrievals.

Gaussian flat line rejection tests
We begin our assessment of the GJ 341b transmission spectrum by quantifying the statistical significance of any potential features in the spectrum relative to the null hypothesis.In this case, the null hypothesis is that the planet possesses a flat, featureless transmission spectrum that is devoid of atmospheric information, as would be the case if the planet does not possess an atmosphere or the atmospheric signal falls below our detection sensitivity (Kreidberg et al. 2014;Lustig-Yaeger & Fu et al. 2023b).
Following the approaches of Moran & Stevenson et al. (2023) and May & MacDonald et al. (2023), we perform two statistical tests on the transmission spectra from each data reduction pipeline.The first test quantifies the (in)consistency of our measurements with the null hypothesis given the observational errors, while the second test quantifies the evidence favoring an agnostic Gaussian absorption feature in the spectrum over the flat null hypothesis.As both tests use slightly different approaches to quantify the significance of one or more features in the spectrum, they are both useful diagnostics.However, these tests are limited to statistical statements about the spectrum as opposed to physical explanations, which require the more detailed modeling presented in the forthcoming subsections.
To understand whether the observed spectrum is consistent with the null hypothesis, we generated 100,000 random spectra on the same wavelength grid subject to the assumption that the underlying true planet spectrum is the null hypothesis and the observational errors are unchanged from those observed.We calculated the reduced chi-squared (χ 2 ν ) for each sample and build up the expected distribution of χ 2 ν under the null hypothesis.That is, if the null hypothesis were true then we would expect that our measured χ 2 ν would be a single sample from that distribution.Thus, we can quantify the probability that the measured χ 2 ν value would be at least as large as the one we observed given the assumption of the null hypothesis (i.e., a p-value).
To understand whether there is evidence of a spectral feature in the observed spectrum, we first fit a flat featureless spectrum to the data and then fit an agnostic Gaussian model to the spectrum.We then compare the Bayesian evidence of the two fits and calculate Bayes factors to estimate which model is preferred, as in Moran & Stevenson et al. (2023) and May & MacDonald et al. (2023).The Gaussian model contains four fitting parameters, which are all constrained by uniform priors.These include the wavelength independent transit depth upon which the Gaussian feature is added ∼U(0.0, 1.0), the wavelength of the Gaussian center ∼U(3.5, 5.5) µm, the width of the feature ∼U(∆x 0 , ∆x 1 ) µm, where ∆x 0 is the minimum spectroscopic bin width and ∆x 1 is the width of the entire wavelength range of the data, and the amplitude of the Gaussian feature ∼U(−10σ y , 10σ y ), where σ y is the standard deviation of the transmission spectrum measurements (i.e., features should not greatly exceed the scatter in the data).
The results of these tests are shown in Table 3.This demonstrates that the Eureka!spectrum is the flattest of the three spectra, ruling out our null hypothesis at only 0.2σ, while the Gaussian feature test favors a flat line over a feature at 2.2σ.The Tiberius and Tswift spectra are less flat.However, in the case of the Tiberius spectrum we believe this is due to the smaller uncertainties in its transmission spectrum (see Section 4.2).For these reductions, the null hypothesis is in weak tension, at 1.3σ for Tiberius and 1.4σ for Tswift.The Gaussian feature test results in 2.3σ evidence for a feature in the Tiberius spectrum (broad and centered near 4.3 µm) and 2.9σ in the Tswift spectrum (width of ∼0.2 µm centered near 4.7 µm).While these look promising at face value, as we see from these tests and demonstrate with our subsequent retrieval analysis (section 4.4), the two candidate features occur at different wavelengths between the Tiberius and Tswift spectra and are therefore less likely to be astrophysical.

Sensitivity to Spectrum Uncertainties
The results of the aforementioned flat line rejection tests depend not just on potential features in the spectrum, but also on the size of the error bars, which are not equivalent for each reduction.The 1σ spectroscopic 100% Photosphere: 3500 K ( 2 =87.1) 56% Photosphere: 3500 K 44% Spot: 3300 K ( 2 =52.9) 54% Photosphere: 3300 K 14% Spot: 3100 K 32% Faculae: 3500 K ( 2 =66.2) GJ 341 Figure 5.A comparison between our measured flux calibrated spectra (black) and superpositions of PHOENIX stellar atmosphere models: photosphere with no active region contribution (purple), photosphere with contribution from spots (orange), and photosphere with contributions from spots and faculae (green).Given the poor χ 2 ν , and lack of contamination evident in our transmission spectra, we do interpret these model results as evidence for stellar contamination in our data.transit depth uncertainties are shown in Figure 7 (lower panel).Since all reductions used approximately the same spectral bin widths, we can readily compare their 1σ uncertainties and apply the errors from one reduction to the others to estimate the role of uncertainties in driving the statistical results of our null hypothesis tests.Table 4 shows χ 2 ν values for flat line fits to each reduction's transmission spectrum data points (rows) with errors interpolated from the other reductions (columns).
Figure 7 (lower panel) and Table 4 show that the differences in uncertainties between reductions leads to sta-tistically significant differences in χ 2 ν .The Tiberius reduction has the smallest uncertainties and when those are applied to the Eureka! and Tswift spectra, those spectra deviate more from a flat line.On the other hand, the Tswift spectrum has χ 2 ν > 1 regardless of which errors are applied, consistent with the larger deviations seen in the Tswift spectrum.The Eureka! spectrum has relatively large errors, but relatively small spectroscopic deviations, which together explain our statistical result and subsequent physical interpretations.
The upper panel of Figure 7 shows the light curve residual RMS and the expected noise for comparison, Note-χ 2 ν is the reduced chi-squared resulting from the best fitting featureless fit (null hypothesis) to the observed spectrum, "pvalue" refers to the probability that χ 2 ν would be at least as extreme as the observed value under the assumption that the null hypothesis is correct and is displayed along with the corresponding "sigma" value, and the "Feature?"column shows the level of confidence in the detection of an agnostic Gaussian absorption feature in the spectrum over the null hypothesis.which includes the photon noise, read noise and dark current.The Tiberius, Tswift, and Eureka!reductions achieve a median precision of 1.03×, 1.05×, and 1.14× the noise limit, respectively.This ordering is different to the white light curve residuals, where Eureka!obtained the lowest residual RMS (Figure 2).The lower RMS of Tiberius and Tswift at the spectroscopic light curve stage comes from those pipelines' use of a common mode correction.Although the Tiberius transit depth uncertainties are lower than the predicted PandExo uncertainties, this is likely due to improved performance relative to the model-based prediction, since the Tiberius reduction still sits above the noise limit (Figure 7, upper panel).

Atmosphere forward models
As with previous targets in program JWST-GO-1981, we compare our transmission spectra of GJ 341b to a set of atmospheric forward models to infer whether the data suggests that the planet possesses an appreciable atmosphere, and if so, its composition (Lustig-Yaeger & Fu et al. 2023b; Moran & Stevenson et al. 2023;May & Mac-Donald et al. 2023).Since GJ 341b has no published mass constraints beyond an upper limit (≤ 3 M ⊕ at 1σ; DiTomasso et al. 2023), we generate forward models using a planetary mass of 0.73 M ⊕ .This mass assumes that the planetary radius of 0.92 R ⊕ is consistent with a rocky bulk composition, as Luque & Pallé (2022) showed that very small (≤ 1 R ⊕ ) planets around M-dwarfs, like GJ 341b, appear to have very tightly constrained densities similar to Earth-like iron-silicate bulk compositions.
Our atmospheric forward models include two cases: 1) CHIMERA thermochemical equilibrium models (Line et al. 2013(Line et al. , 2014) ) from 100-1000× solar metallicity, using a parameterized temperature-pressure profile (Guillot 2010) with an equilibrium temperature of 540 K and opacity contributions from H 2 O, CH 4 , CO, CO 2 , NH 3 , HCN, H 2 S, H 2 , and He; and 2) PICASO (Batalha et al. 2019) models of single-gas atmospheres with an isothermal temperature-pressure profile at 540 K.The PICASO models include atmospheres of pure 1 bar H 2 O, CO 2 , or CH 4 , as well as a flat line model representative of no atmosphere or a uniform aerosol layer.For all model configurations, we compute synthetic transmission spectra with PICASO's radiative transfer routine using molecular opacities resampled from R = 60,000 (Batalha et al. 2020), rebinned to the resolution of each data reduction.We compute best fits by calculating the χ 2 ν between each reduction and each forward model, with 57, 54, or 55 degrees of freedom (dof) for the Tiberius, Tswift, and Eureka!reductions, respectively.We include an offset in relative transit depth to the mean of each reduction in our fitting.Our results for Eureka! are shown in Figure 8, while those for Tiberius and Tswift are shown in Figures 14 and 15.We include the complete set of statistical constraints for each reduction compared to each forward model case in Table 5.
Across all three reductions, hydrogen-dominated or methane atmospheric models do not fit the data.For the Tswift reduction, hydrogen-dominated atmospheres as metal-rich as 1000× solar can be ruled out to 3σ, while this constraint is reduced to 350× solar for Tiberius and Eureka!.Such hydrogen-rich atmospheres would be unexpected to persist against escape, so this result is unsurprising.Methane atmospheres of 1 bar are disfavored at 3σ, 2.5σ, or 1σ between Tswift, Tiberius, and Eureka!.
On the other hand, water atmospheres or a flat line -which could stem either from an airless world or from a hazy world -are consistent with each reduction, with χ 2 ν ≲ 1.Though water is a ubiquitous molecule and common in warm planetary atmospheres, we are unable to rule it out here since water has no strong absorption bandheads and instead acts as a continuum-like absorber from 4-5 µm.Both Eureka! and Tiberius also are consistent (χ 2 ν ≲ 1) with a carbon dioxide atmosphere, while Tswift disfavors this scenario to 2σ.
Despite the excellent agreement between each data reduction (Figure 6), there are obvious differences in goodness of fit between atmospheric models.This occurs due to 1) the planet's small size and temperature, which results in a very small atmospheric scale height, and 2) the ratio of the sub-Earth radius (R P =0.92 R ⊕ ) to the relatively large radius of the star for an Mdwarf (R * =0.506 R ⊙ ), resulting in a very shallow transit.Given these physical parameters, our 350× solar, hydrogen-dominated atmosphere model has a peak feature size of 40 ppm while a carbon dioxide atmosphere model has a peak feature size of only 15 ppm in the NIRCam/F444W bandpass.Even with the extremely high precision obtained by our program (average uncertainties of 15 ppm per 0.02 µm bin), such small feature sizes mean that minute differences between reductions driven by random noise or data treatments drive divergent atmospheric interpretations.Overall, either no atmosphere, a hazy atmosphere, or an atmosphere containing a species that does not have prominent molecular bands across the NIRCam/F444W bandpass -e.g., a water-dominated atmosphere -provide the best explanation of the data across the three reductions compared to our forward models.

Atmosphere retrievals
To examine a broader range of atmospheric parameter space that is (in)consistent with our measurements, we ran a set of retrievals using the Spectral Mapping Atmospheric Radiative Transfer Exoplanet Retrieval code (smarter; Lustig-Yaeger et al. 2022, 2023a;Lustig-Yaeger & Fu et al. 2023b).We use the same smarter modeling approach and assumptions as described in Lustig-Yaeger & Fu et al. (2023b), except slightly modified to account for the different planet parameters.We run separate retrievals on each of the final coadded spectra from the Eureka!, Tiberius, and Tswift reductions.
Our nominal retrieval setup for GJ 341b uses 11 free parameters.We include the following six spectroscopically active gases that have (or partially have) absorption bands in the NIRCam F444W wavelength range and are plausible candidate gases for secondary atmospheres at temperatures at or near GJ 341b's equilibrium temperature of 540 K: H 2 O, CH 4 , CO, CO 2 , N 2 O, and O 3 .Each gas is fitted in terms of the log 10 evenly-mixed volume mixing ratio (VMR) with a flat prior, VMR ∼ U(−12, 0).The bulk atmospheric constituent is left agnostic by fitting for the mean molecular weight of the primary constituent using a flat prior, µ ∼ U(2.5, 50) g/mol, which ranges from a low-mass mix of H 2 +He to heavier than CO 2 (44 g/mol) and O 3 (48 g/mol).Although aerosols are not explicitly included in our retrievals, we fit for the opaque planet radius using a uniform prior R p ∼ U(0.83, 1.01)R ⊕ , which could be due to a cloud deck, the solid planetary surface (Benneke & Seager 2012), additional continuum opacity, or the critical refraction limit (Bétrémieux & Kaltenegger 2013;Misra et al. 2014).We fit for the pressure at the opaque surface using a uniform prior, P 0 ∼ U(−1, 6) Pa.An isothermal temperature structure is assumed  with a uniform prior, T ∼ U(300, 700) K. Finally, we allow the planet mass to vary uniformly between M p ∼ U(0.3, 2.2)M ⊕ .
Our smarter retrievals use the dynesty (Speagle 2020) nested sampling algorithm (Skilling 2004) with 600 live points.We run the model until the contribution to the total evidence from the remaining prior volume drops below dlogz=0.075.
Figure 9 shows a summary of our retrieval results, comparing top-level findings for each data reduction pipeline.Figures 16,17,and 18 (in Appendix E) show results for the Eureka!, Tiberius, and Tswift pipelines, respectively.Based on the χ 2 ν values from the median retrieved spectrum, the retrievals are over-fitting the spectrum, and therefore likely fitting noise.Nonetheless, they reveal important aspects of our observational sensitivity to the possible atmosphere of GJ 341b as well as the agreement/disagreement between the reduction pipelines.As previously discussed, the Eureka!reduction appears the most flat and as a result the retrieval rules out a limited slice of parameter space that corresponds to only the thickest (reference pressures P 0 > 10 4 Pa) and most extended atmospheres, with high abundances of absorbing gases, that we explored.
On the other hand, the Tiberius and Tswift reductions tell a different story due to their low significance detection of a spectroscopic feature (Table 3), but they do not agree on the wavelength/molecule responsible for the feature, implying that these are not real astrophysical signals.The retrieval results for Tiberius favor high abundances of CO 2 (−2.1 +1.2 −1.9 dex) with moderate abundances of N 2 O (−4.3 +2.1 −3.9 dex) and low abundances of O 3 (−8.3+2.5  −2.7 dex; 2σ upper limit of −3.6 dex).The retrieval results for Tswift show the opposite, favoring low abundances of CO 2 (−8.7 +2.6  −2.1 dex; 2σ upper limit of −3.6 dex) and N 2 O (−8.2 +2.8 −2.4 dex; 2σ upper limit of −2.7 dex) and high abundances of O 3 (−1.8+1.0 −1.8 dex).These O 3 and CO 2 constraints differ at ∼1.5σ.Our Bayesian model comparisons (using the Bayesian evidence ratios computed by dynesty) find that the Tiberius reduction implies a weak detection of CO 2 (3.1σ), but no O 3 , while the Tswift reduction instead leads to a weak detection of O 3 (2.9σ),but no CO 2 .Conversely, the Eureka!reduction favors no detected gases at >1σ significance.However, both Tiberius and Tswift retrievals levy similar constraints on the reference pressure, isothermal temperature, reference radius, planet mass, and mean molecular weight which are all in agreement at better than 1σ.These findings suggest that different noise instances with otherwise similar noise properties dominate the character of each reduced spectrum, as opposed to planetary characteristics.Furthermore, the statistical significance of results from any one pipeline must be interpreted as lower than the assessment of that spectrum in isolation would suggest.

DISCUSSION AND CONCLUSIONS
We have measured the 3.9-5.0µm transmission spectrum of the sub-Earth sized exoplanet GJ 341b using JWST's NIRCam instrument.Our transmission spectrum is the result of co-adding the transmission spectra measured from three separate transit observations and represents the first observation of this planet's atmosphere.We reduced the data using three independent pipelines (Eureka!, Tiberius and Tswift) and performed several tests to determine the significance of an atmosphere detection.

The NIRCam/F444W data quality
Our NIRCam/F444W data is dominated by an exponential ramp that varies between the three visits (Figure 1).Despite fitting for this ramp, there is residual red noise in the white light curve data (Figure 2).This red noise becomes apparent at timescales between ∼ 1 − 30 minutes.We attempted to remove this red noise by using up to ten engineering telemetry timeseries ('mnemonics') in our systematics model (Figure 3).Despite encouraging signs, we were unable to significantly improve the residual red noise.We believe that the JWST engineering mnemonics do present an opportunity to detrend JWST spectrophotometry in general, however, we leave further attempts to do so to future work.Fortunately, this red noise is not significant in the spectroscopic light curves with our light curves reaching between 1.03 and 1.14× the noise limit (Fig. 7).As a result, we were able to obtain a very precise transmission spectrum from our three coadded visits, with an average uncertainty of 15 ppm in 0.02 µm bins.For these reasons, we recommend the use of NIRCam/F444W for studies of bright exoplanet systems.While NIRCam also provides us with short wavelength photometry (as used successfully in the ERS program, Ahrer et al. 2023), our extraction of this here shows that is not always beneficial.

GJ 341b's transmission spectrum
The main conclusion we draw from our high precision transmission spectrum of GJ 341b is that the planet has not retained a low mean molecular weight atmosphere.This is consistent with expectations from the cosmic shoreline (Zahnle & Catling 2017) that small planets like GJ 341b should be stripped of any primordial H/Hedominated atmosphere.Digging deeper, we find that the very high precision of our data, combined with the small expected feature amplitudes, leads our statistical tests, forward models and retrievals to slightly different conclusions.
In our null hypothesis tests, the Eureka!spectrum was consistent with being flat to within 0.2σ while the Tiberius and Tswift spectra were less consistent at 1.3σ and 1.4σ, respectively.However, we believe that the relative consistency with being flat is driven by the relative size of the uncertainties for each transmission spectrum.Repeating the null hypothesis tests but switching the uncertainties between the reductions, we find that Tiberius is inherently the flattest spectrum, followed by Eureka! and then Tswift.
Our Gaussian feature tests, using the default uncertainties for each reduction, find 2.2σ evidence against a feature for Eureka! and 2.3σ and 2.9σ evidence for features in the Tiberius and Tswift spectra, respec-tively.However, our atmospheric retrievals find that different species are responsible for these candidate features: CO 2 in the Tiberius spectrum at 3.1σ and O 3 in the Tswift spectrum at 2.9σ, with a relative difference between the CO 2 and O 3 abundance constraints of ∼1.5σ.These differences suggest that these 'features' are not real astrophysical signals.
Our forward model analysis finds that each reduction's spectrum rejects a < 285× solar metallicity atmosphere to ≥ 3σ, with the Tswift spectrum rejecting a < 1000× solar metallicity atmosphere.Instead, the forward models find our transmission spectra are consistent with no atmosphere, a hazy atmosphere, or an atmosphere containing a species that does not have prominent molecular bands across the NIRCam/F444W bandpass, e.g.H 2 O. Furthermore, CH 4 -dominated atmospheres are ruled out to 1σ-3σ depending on the reduction, while CO 2 atmospheres are consistent with the Eureka! and Tiberius spectra to < 1σ and inconsistent with the Tswift spectrum to > 2σ.
At the radius and equilibrium temperature of GJ 341b, atmospheres dominated by CO 2 , H 2 O or O 2 /O 3 , resulting from post-runaway atmospheres, are all plausible and even predicted (Luger & Barnes 2015;Gao et al. 2015;Schaefer et al. 2016;Lincowski et al. 2018;Kite & Schaefer 2021;Krissansen-Totton & Fortney 2022).Therefore, the single gas end member atmospheres that result from our forward modelling, and that are consistent with our retrievals, are in line with expectations.However, since a rocky exoplanet atmosphere has yet to be detected and characterized, the possibilities are largely restricted to theoretical predictions and the range of possible temperatures, pressures and gas abundance combinations that are consistent with the observations is vast.Furthermore the fact that each pipeline prefers a different gas demonstrates that the evidence for these gases is not due to real astrophysical signals.
Our study of this shallow transit (240 ppm) demonstrates that small differences between pipelines' spectra can lead to differences in the conclusions drawn from statistical tests, atmospheric forward models and retrievals.Therefore, it is imperative that such studies use multiple independent pipelines, two or more visits, and several statistical tests to determine the robustness of future claims of atmospheric absorption for small exoplanets.

The Cosmic Shoreline for M-dwarfs
GJ 341b is the fourth planet to be studied as part of our reconnaissance program (PID:1981) to test the existence of an M-dwarf Cosmic Shoreline (Zahnle & Catling 2017).Our results so far are consistent with expectations (Figure 10).GJ 341b (this work) and LHS 475b (Lustig-Yaeger & Fu et al. 2023b) sit on the dry/airless side of the shoreline, which is consistent with our transmission spectra.GJ 486b (Moran & Stevenson et al. 2023) and GJ 1132b (May & MacDonald et al. 2023) sit at the edge of the shoreline, where atmospheres and atmosphere-less are both possible outcomes.GJ 486b shows water in its transmission spectrum, although this could be stellar contamination.GJ 1132b paints an inconsistent picture between the two visits' transits, making inferences more challenging.Our remaining target (TRAPPIST-1h) is likely to have an atmosphere based on expectations from the shoreline.Our upcoming transit observations will shed light on the existence of its atmosphere.Taken together, our program has demonstrated the tight constraints that can be placed on Earth-sized exoplanet atmospheres.
Looking beyond our program to other JWST studies of Earth-sized exoplanets around M-dwarfs, TRAPPIST-1b and c are both consistent with having no atmosphere (Greene et al. 2023;Lim et al. 2023;Zieba et al. 2023).TRAPPIST-1b is close to the Cos-mic Shoreline and thus the lack of an atmosphere is not in tension with expectations.TRAPPIST-1c, on the other hand, sits on the atmosphere side.However, the shoreline presented in Figure 10 is plotted in terms of the relative instellation.The XUV cosmic shoreline would sit below this (Zahnle & Catling 2017) and may be consistent with the TRAPPIST-1c result.If the results of no atmospheres for these planets hold, then this could be emerging evidence that XUV dominates the cosmic shoreline in these late M-dwarf systems, consistent with past escape predictions (Wheatley et al. 2017;Dong et al. 2018).Only by combining results from wider and deeper surveys, including ongoing JWST programs (e.g., PID: 2512), will we be able to confidently determine the boundaries of a cosmic shoreline for M-dwarfs, if indeed one exists.

Figure 6 .
Figure6.GJ 341b's transmission spectrum as measured by Tiberius (blue), Eureka!(orange) and Tswift (green).These spectra are the the weighted means of the three visits' spectra.The Eureka! spectrum has been offset by -15 ppm to make the median transit depths equal to Tiberius and Tswift between 4.0 and 4.8 µm.

Figure 7 .
Figure 7.Light curve residual RMS (upper panel) and transit depth uncertainties (lower panel) for each reduction (color lines).The expected noise (photon + readnoise + dark current) calculated from the data is shown in the upper panel (black line) and the median RMS relative to the expected noise is listed in the legend for each reduction.The predicted noise from PandExo (Batalha et al. 2017) is shown in the lower panel (black line).

Figure 8 .
Figure8.Atmospheric forward models compared to the Eureka!reduction from JWST/NIRCam.Because of the planet's small size and low equilibrium temperature, only hydrogen-dominated or methane atmospheres are strongly ruled out or disfavored by the data.Atmospheres of water, carbon dioxide, or an opaque haze layer -or no atmosphere -are all statistically consistent with the data to within ∼1σ across all three reductions.

Figure 9 .
Figure 9.Comparison of retrieval results on data reduced by the Eureka!, Tiberius, and Tswift pipelines, respectively (rows of subplots).For each reduction's retrieval, the median spectral fit is shown in the left panel, the 1D marginalized CO2 posterior is shown in the middle panel, and the 1D marginalized O3 posterior is shown in the right panel.Although either CO2 or O3 are weakly detected in retrievals of the Tiberius and Tswift spectra, respectively, (the Eureka!result is unconstrained) these findings are incompatible with one another.Thus, relatively small differences in the reductions can cause comparatively large changes in the subsequent physical interpretation, including spurious molecular detections.

Facilities:
JWST(NIRCam) All of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute.The specific observations analyzed can be accessed via DOI.Software: Astropy (Astropy Collaboration et al. 2013, 2018), batman (Kreidberg 2015), CHIMERA (Line & Yung 2013; Line et al. 2014), Dynesty (Speagle 2020), emcee (Foreman-Mackey et al. 2013), Eureka!(Bell et al. 2022), ExoCTK (Bourque et al. 2021), ExoTiC-LD (Grant & Wakeford 2022) FIREFLy (Rustamkulov et al. 2022), Forecaster (Chen & Kipping 2017), IPython c SLC corresponds to the spectroscopic light curves.A. THE REDUCTION AND FITTING APPROACHES TAKEN BY Eureka!, Tiberius AND Tswift B. THE LIGHT CURVES AND RESIDUALS FROM EACH PIPELINE C. TRANSMISSION SPECTRA MEASURED FOR EACH VISIT The transmission spectra measured for each visit by each reduction are shown in Figure 13.

Figure 13 .
Figure 11.The light curves and best fit models resulting from the Tiberius pipeline.The columns show, from left to right: transit 1, transit 2, transit 3. The rows show, from top to bottom: 1) the white light curve and residuals, 2) the spectroscopic light curves, 3) the best-fit models to the spectroscopic light curves, 4) the residuals from the spectroscopic light curve fitting.The vertical band at the left hand edge of the panels in rows 2-4 correspond to the 600 integrations that were masked during the spectroscopic light curve fitting stage.

Figure 14 .
Figure 14.Atmospheric forward models compared to the Tiberius reduction from JWST/NIRCam.
Figure3.The improvement in residual red noise in the white light curves by adding an additional mnemonic / telemetry timeseries to the systematics model.The blue, orange and green lines indicate the residuals from visit 1, 2 and 3's white light curves.The fiducial model is a linear-intime polynomial.Each x-tick label indicates an additional mnemonic is added.The dashed gray line indicates the white noise level that the residual noise would achieve if no red noise is present.

Table 2 .
The ten engineering mnemonics that led to the largest improvement in the cumulative residual red noise from all three white light curve fits.These are ordered from largest improvement (top) to least improvement (bottom).With regards to the mnemonics, 'INRC' is NIRCam telemetry reported by ISIM, the Integrated Science Instrument Module, 'IGDP' are 'ground data points' computed from downlinked engineering telemetry, and 'TEMP' corresponds to temperature.We understand that 'A' and 'B' correspond to module A and module B, 'ACE' corresponds to Application Specific Integrated Circuit (ASIC) control electronics, and 'LV' corresponds to low voltage.We define the residual red noise in the text.

Table 3 .
Is it Flat?

Table 4 .
Reduced-χ 2 Values for a flat line fit when exchanging errors between reductions Tiberius Errors Eureka!Errors Tswift Errors

Table 6 .
Comparison of the different reduction and analysis choices used by the three pipelines, along with the resulting noise properties.The goal of group-level background subtraction is to remove 1/f noise.bWLC corresponds to the white light curves. a