The search for rocky planet atmospheres with JWST has focused on planets transiting M dwarfs. Such planets have favorable planettostar size ratios enhancing the amplitude of atmospheric features. Since the expected signal strength of atmospheric features is similar to the singletransit performance of JWST multiple observations are required to confirm any detection. Here we present two transit observations of the rocky planet GJ 1132 b with JWST NIRSpec G395H covering 2.85.2 m. Previous Hubble Space Telescope WFC3 observations of GJ 1132 b were inconclusive with evidence reported for either an atmosphere or a featureless spectrum based on analyses of the same data set. Our JWST data exhibit substantial differences between the two visits. One transit is consistent with either an HSUB2SUBOdominated atmosphere containing 1 CHSUB4SUB and trace NSUB2SUBO  chi nu 21.13  or stellar contamination from unocculted starspots  chi nu 21.36 . However the second transit is consistent with a featureless spectrum. Neither visit is consistent with a previous report of HCN. Atmospheric variability is unlikely to explain the scale of the observed differences between the visits. Similarly our outoftransit stellar spectra show no evidence of changing stellar inhomogeneity between the two visitsobserved 8 days apart only 6.5 of the stellar rotation rate. We further find no evidence of differing instrumental systematic effects between visits. The most plausible explanation is an unlucky random noise draw leading to two significantly discrepant transmission spectra. Our results highlight the importance of multivisit repeatability with JWST prior to claiming atmospheric detections for these small enigmatic planets.
['Exoplanet atmospheres', 'Exoplanet atmospheric composition', 'Transmission spectroscopy', 'Astronomy data reduction', 'Planetary atmospheres', 'Stellar atmospheres', 'Starspots', 'Super Earths', 'Infrared spectroscopy', 'Extrasolar rocky planets', '487', '2021', '2133', '1861', '1244', '1584', '1572', '1655', '2285', '511', 'Astrophysics - Earth and Planetary Astrophysics']
Double Trouble Two Transits of the SuperEarth GJ 1132 b Observed with JWST NIRSpec G395H
{'Double Trouble: Two Transits of the Super-Earth GJ 1132 b Observed with JWST NIRSpec G395H': '<!-- image -->  \nKristin S. Sotzen 1, 2  \n1 Johns Hopkins APL, Laurel, MD 20723, USA  \n2 Consortium on Habitability and Atmospheres of M-dwarf Planets (CHAMPs), Laurel, MD, USA  \n3 Department of Astronomy, University of Michigan, Ann Arbor, MI 48109, USA  \n4 NHFP Sagan Fellow  \n5 Department of Earth & Planetary Sciences, Johns Hopkins University, Baltimore, MD, USA  \n6 Department of Planetary Sciences and Lunar and Planetary Laboratory, University of Arizona, Tuscon, AZ, USA  \n7 School of Physics, HH Wills Physics Laboratory, University of Bristol, Bristol, UK  \n8 University of Maryland, Baltimore County, MD 21250, USA  \n9 NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA  \n10 Department of Physics & Astronomy, Johns Hopkins University, Baltimore, MD, USA  \n11 NASA Ames Research Center, Moffett Field, CA 94035, USA  \n12 Department of Physics, Imperial College London, Prince Consort Road, London, SW7 2AZ, UK  \n13 Center for Astrophysics | Harvard & Smithsonian, 60 Garden St, Cambridge, MA 02138, USA  \n14 Space Telescope Science Institute, Baltimore, MD 21218, USA  \n15 Carnegie Earth & Planets Laboratory, Washington, DC, 20015 USA  \n16 Department of Physics, Utah Valley University, Orem, UT, 84058 USA', 'ABSTRACT': 'The search for rocky planet atmospheres with JWST has focused on planets transiting M dwarfs. Such planets have favorable planet-to-star size ratios, enhancing the amplitude of atmospheric features. Since the expected signal strength of atmospheric features is similar to the single-transit performance of JWST , multiple observations are required to confirm any detection. Here, we present two transit observations of the rocky planet GJ 1132 b with JWST NIRSpec G395H, covering 2.8-5.2 µ m. Previous HST WFC3 observations of GJ 1132 b were inconclusive, with evidence reported for either an atmosphere or a featureless spectrum based on analyses of the same dataset. Our JWST data exhibit substantial differences between the two visits. One transit is consistent with either a H 2 O-dominated atmosphere containing ∼ 1% CH 4 and trace N 2 O ( χ 2 ν = 1 . 13) or stellar contamination from unocculted starspots ( χ 2 ν = 1 . 36). However, the second transit is consistent with a featureless spectrum. Neither visit is consistent with a previous report of HCN. Atmospheric variability is unlikely to explain the scale of the observed differences between the visits. Similarly, our out-of-transit stellar spectra show no evidence of changing stellar inhomogeneity between the two visits - observed 8 days apart, only 6.5% of the stellar rotation rate. We further find no evidence of differing instrumental systematic effects between visits. The most plausible explanation is an unlucky random noise draw leading to two significantly discrepant transmission spectra. Our results highlight the importance of multi-visit repeatability with JWST prior to claiming atmospheric detections for these small, enigmatic planets.', '1. INTRODUCTION': "The quest to detect atmospheres on rocky exoplanets requires pushing our observatories to their limits. Even for rocky planets transiting M dwarfs - the most promising for atmospheric detections - transmission spectra features have an expected amplitude of ≲ 20 ppm. Atmospheric features of such planets are thus perilously close to the pre-launch expected noise floor of JWST instruments ( ∼ 20 ppm for NIRISS, ∼ 9 ppm for NIRCam, and < 10 ppm for NIRSpec; Greene et al. 2016; Schlawin et al. 2021; Rustamkulov et al. 2022). This comparable size of features and noise therefore calls for extra care before claiming an atmospheric detection. In particular, the repeatability of any signal between visits with the same instrument is critical to confirm that observed features are astrophysical in nature.  \nThere has been no definitive or non-disputed detection of an atmosphere on a rocky exoplanet to date. Rocky bodies in the solar system exhibit a 'cosmic shoreline', which divides bodies with and without atmospheres according to the prevalence of atmospheric escape processes (Zahnle & Catling 2017). Such a classification scheme may be a logical starting point for selecting promising rocky exoplanets for atmospheric characterization. However, planets orbiting M dwarfs may exhibit a significantly different cosmic shoreline, or no cosmic shoreline at all, due to the higher stellar activity and extreme-UV (1-912 ˚ A) flux levels compared to the Sun, which can strip away planetary atmospheres (e.g., Owen & Jackson 2012; Becker et al. 2020; Dong et al. 2020; do Amaral et al. 2022). In fact, recent secondary eclipse observations with the Mid-Infrared Instrument (MIRI) on JWST of TRAPPIST-1 b and TRAPPIST-1 c, two rocky M-dwarf planets, are consistent with no (or minimal) atmosphere (Greene et al. 2023; Zieba et al. 2023).  \nSeveral JWST Cycle 1 programs are surveying rocky planets orbiting M dwarfs, aiming to assess the survivability of their atmospheres (e.g., JWST GO #1981, 2512, 2589). Through JWST GO #1981 (PIs: Stevenson & Lustig-Yaeger), we are observing five rocky planets around M dwarfs with orbital and planetary properties close to the Solar System's cosmic shoreline. Previous results from this program include two transits of LHS 475 b, which are consistent with no atmosphere or a high-altitude cloud deck (Lustig-Yaeger & Fu et al. 2023), and two transits of GJ 486 b, which are consistent with a H 2 O-dominated atmosphere or a spectrum contaminated by stellar activity (Moran & Stevenson et  \nal. 2023). Here we present transmission spectra observations for the third planet in the program, GJ 1132 b.  \nGJ 1132 b (Berta-Thompson et al. 2015) is a rocky super-Earth (1.30 R ⊕ , 1.66 M ⊕ , T eq = 529K; Bonfils et al. 2018) orbiting an M dwarf (0.2105 R ⊙ , 3261 K; Bonfils et al. 2018). With a favorable planet-star radius ratio, corresponding to a transit depth of ∼ 0.3%, GJ 1132 b has previously been suggested as a good target for atmospheric characterization (e.g. Schaefer et al. 2016). Assuming a representative atmospheric mean molecular mass of µ = 10AMU, GJ 1132 b should have transmission spectra features of ∼ 20 ppm.  \nPrevious transmission spectra observations of GJ 1132 b have a colorful history of claimed atmospheric detections. Southworth et al. (2017) observed nine transits with the MPG 2.2 m telescope and suggested that deeper transit depths in the z and K bands were caused by a H 2 -dominated atmosphere with H 2 O and/or CH 4 . Diamond-Lowe et al. (2018) revisited GJ 1132 b with five optical transits (0.64-1.04 µ m) with the Magellan Clay telescope LDSS-3C instrument, finding a featureless transmission spectrum ruling out the previously-claimed atmosphere. Swain et al. (2021) analyzed five near-infrared (1.1-1.7 µ m) Hubble Space Telescope Wide Field Camera 3 (WFC3) transits ( HST GO #14758, PI: Berta-Thompson), finding evidence for a spectral slope and a feature near 1.53 µ m suggestive of a H 2 -dominated atmosphere with aerosols, HCN, and CH 4 . However, using the same HST data, Mugnai et al. (2021) and Libby-Roberts et al. (2022) do not find these features and instead prefer featureless near-infrared spectra.  \nIn this Letter, we present JWST transmission spectra observations of GJ 1132 b. Our analysis provides a cautionary tale for the challenge of confirming astrophysical signals when searching for rocky exoplanet atmospheres. In Section 2, we describe the observations. In Section 3, we overview the three independent analysis pipelines used to extract the transmission spectrum. Section 4 describes our interpretation. Finally, we discuss the implications of these results in Section 5.", '2. JWST OBSERVATIONS OF GJ 1132 b': 'We observed two transits of GJ 1132 b with the JWST Near Infrared Spectrograph (NIRSpec, Jakobsen et al. 2022; Birkmann et al. 2022) G395H instrument on 2023 February 25 and 2023 March 5 as a part of GO #1981 (PIs: Stevenson & Lustig-Yaeger). Our data have a spectral resolving power of λ/ ∆ λ ∼ 2,700 from 2.85.2 µ m. Each observation lasted 3.06 hrs, resulting in 814 integrations each with 14 groups up the ramp. The  \nobservations were designed to maximize observing efficiency while remaining below the ∼ 80% full well threshold to avoid the worst impacts of detector non-linearity.', '3. DATA REDUCTION': 'To ensure reproducibility of our results, we perform three independent data analyses with the Eureka! (Bell et al. 2022), FIREFLy (Rustamkulov et al. 2022, 2023), and ExoTiC-JEDI (Alderson et al. 2022, 2023) pipelines. Previous analyses within JWST GO #1981 have shown Eureka! and FIREFLy to agree well for small planetary signals (Lustig-Yaeger & Fu et al. 2023; Moran & Stevenson et al. 2023). While Eureka! and ExoTiC-JEDI have been shown to agree well for larger planetary signals (e.g. Alderson et al. 2023), here we compare them for signals pushing detection limits. Below, we provide a high-level overview of the data reduction steps from each pipeline. Figure 1 shows the final R ∼ 100 spectra from all three pipelines for both visits.', '3.1. Eureka!': 'The Eureka! package reduces JWST time-series data starting from uncal JWST data through light curve fitting. Stages 1 and 2 of the Eureka! pipeline are primarily a wrapper for Stages 1 and 2 of the jwst pipeline (Bushouse et al. 2022), while also implementing several custom steps, most importantly custom group-level background subtraction. This removes the striping due to 1/f noise at the group level to significantly improve precision (e.g., Rustamkulov et al. 2022; Lustig-Yaeger & Fu et al. 2023). This step first identifies the center of light in each column and masks values within 8 pixels on either side. We estimate the background as the mean of all remaining pixels in a column with a 3 σ outlier threshold. We skip the jump step detection in Stage 1, as we find it only adds noise to the extracted light curves. In Stage 2, we skip the flat field and photom steps, as they result in a conversion to physical flux units that is unnecessary for the relative flux measurements we require. These steps also add noise to the resulting light curves due to the current level of accuracy provided by the available detector flat fields. Additionally, the limited region of the detector that is converted in these steps also worsens the precision on background removal.  \nStage 3 of the Eureka! pipeline performs spectral extraction after a second round of background subtraction. We use optimal spectral extraction (Horne 1986) after correcting the curvature of the trace by measuring the center of light in each column and rolling each column by an integer pixel value to align the trace along the same pixel. The second round of background subtraction considers only the region more than 10 pixels from  \nthe trace, with an outlier threshold of 3 σ in the spatial direction and 10 σ in the temporal direction. When constructing the median frame for optimal spectral extraction, we employ an outlier rejection threshold of 10 σ for NRS1 and 15 σ for NRS2. For spectral extraction, we use an aperture half-width size of 2 pixels on either side of the center pixel (for a total of 5 pixels), which is selected to achieve the best precision possible by minimizing background noise while maximizing the stellar flux extracted. During spectral extraction, we use an outlier rejection threshold of 7 σ for NRS1 and 19 σ for NRS2. In Stage 4, we generate light curves at the native pixel resolution and at a lower resolution of R ∼ 100.  \nTo determine the best orbital parameters in Stage 5 of Eureka! , we first perform joint white light curve fitting across both visits but independently for the two detectors (i.e., NRS1 Visits 1 and 2 are fit jointly). We fit for a linear trend in time combined with a transit function ( batman ; Kreidberg 2015). Limb darkening is fixed to quadratic values obtained with the ExoTiC-LD package (Grant & Wakeford 2022) using 3D stellar models (Magic et al. 2015) and assuming stellar values of T eff = 3261K, log g = 5.02 (both from Stassun et al. 2019), and a metallicity [Fe/H] = -0.12 (Berta-Thompson et al. 2015). We fit for R p /R ⋆ , the center of transit, a/R ⋆ , orbital inclination, and the linear term of the temporal ramp. The planet orbital period is fixed to 1.628931 days (Bonfils et al. 2018). We then adopt the weighted mean of the fitted orbital parameters from the NRS1 and NRS2 joint visit fits as our orbital parameter solution. The resulting best fit values are given in Appendix A.3, Table 2, and are held constant in all spectroscopic light curve fits. The only free parameters are then R p /R ⋆ and the linear temporal ramp terms for both the native pixel and R ∼ 100 resolutions. All fits use emcee (Foreman-Mackey et al. 2013) and are run sufficiently long to ensure chain convergence. 1', '3.2. FIREFLy': "The FIREFLy package undertakes the complete process of JWST time-series data analysis from uncal data files through spectroscopic light curve fitting. We use Stages 1 and 2 of the jwst reduction pipeline for grouplevel and integration-level detector and instrument corrections, respectively. In Stage 1, the data quality initialization and saturation steps are first applied to the uncal fits files. Unlike the reduction of GJ 486 b (Moran & Stevenson et al. 2023), we do not apply a custom superbias scaling step. In this previous study, we inves-  \nFigure 1. Transmission spectra of GJ 1132 b from all three data reduction pipelines and both visits. Upper panel: the Eureka! (yellow), FIREFLy (pink), and ExoTiC-JEDI (blue) reductions are shown for both visits at R ∼ 100. Also shown are the differences between reductions (solid, dashed, and dotted gray lines) in units of σ , demonstrating the agreement between pipelines within a single visit. A ± 1 σ shaded region is overlaid. Bottom panel: the Eureka! reduction for both visits, with the NIRSpec detector gap denoted by the shaded gray region. The shaded yellow regions show two important wavelength ranges that differ the most between the visits (see Sections 4.1 and 4.4.1 for further discussion on these differences).  \n<!-- image -->  \ntigated how the bias levels change over the course of the observation. We found that subtracting a scaled superbias file at the group level improved the agreement between the white light curve depths between the two detectors (i.e., helped account for a possible offset between detectors). The scaling factor was calculated by finding the median background per column in each group, dividing this by the median per column superbias value, and averaging all columns to get a single scaling factor per group. For GJ1132b, we instead elect to use the single default jwst superbias file, as this results in the most consistent white light curve transit depths between the four datasets (NRS1 and NRS2 for Visits 1 and 2). We do, however, implement a custom background subtraction to reduce 1/f noise at the group-level in Stage 1. We then apply the reference pixel correction and linearity step while skipping the dark current step. We also skip the jump step, which is only applied in the FIREFLy pipeline if the number of groups per integration is larger  \nthan 25, as fewer groups per integration (here, 14) lowers the risk of cosmic ray hits. After ramp fitting and the gain-scale step, in Stage 2 we use only the assign WCS step (skipping the flat field step) before proceeding to FIREFLy 's stellar extraction.  \nFor stellar extraction, we first clean bad pixels using lacosmic (van Dokkum 2001). We determined the bad pixel map by flagging pixels with sharp variance spikes and manually checking known bad pixels in NIRSpec G395H. We apply another background subtraction at the integration level, measure the x- and y-shifts, and finally extract the 1D stellar spectrum. We use a precalculated trace and aperture full-width of 5.93 pixels (optimized from previous NIRSpec G395H observations; e.g., Lustig-Yaeger & Fu et al. 2023, Moran & Stevenson et al. 2023) and compute a box extraction.  \nWe fit GJ 1132 b's light curves with batman (Kreidberg 2015), both at the native pixel resolution and at R ∼ 100. For the R ∼ 100 case, we trim the first 150  \ncolumns of NRS1. To fit the white light curve, we fix the orbital period to 1.629 days (Bonfils et al. 2018) and fit for R p /R ⋆ , the mid-transit time T 0 , a/R ⋆ , the impact parameter b , the quadratic limb darkening coefficients, and several systematics parameters. We fit each of the four datasets individually, using the Bayesian Information Criterion (BIC) to determine the best-fit systematics model in each case. Therefore, different models can be used for different datasets. Specifically, we use a linear ramp in time for NRS1 (both visits) and NRS2 Visit 1. For NRS2 Visit 2, we use only the y-shift. After light curve fitting using these optimized systematics, we refit all four datasets using a weighted average of a/R ⋆ , b , and the limb darkening coefficients.  \nTo extract the transmission spectrum, we fit the spectroscopic light curves (holding the orbital parameters and limb darkening coefficients fixed to their white light curve values) with only R P /R ⋆ and the corresponding systematics as free parameters. While FIREFLy utilizes emcee (Foreman-Mackey et al. 2013) when fitting the white light curve, the spectroscopic fits use least-squares fitting with lmfit (Newville et al. 2014). This method is much faster and does not impact the resulting spectrum. We extensively compared both techniques at the native pixel level and at R ∼ 100 and found no meaningful changes when using least-squares fitting.", '3.3. ExoTiC-JEDI': 'The Exoplanet Timeseries Characterisation - JWST Extraction and Diagnostic Investigator ( ExoTiC-JEDI ) package 2 performs a full extraction, reduction, and analysis of JWST time-series data from JWST uncal files to light curve fitting. For NIRSpec G395H observations, we treat the NRS1 and NRS2 datasets independently. We use the jwst pipeline tools to perform linearity, dark current, and saturation corrections, with the jump detection threshold set to 15, and a custom destriping routine to remove 1/f noise at the group level using a second-order polynomial and a threshold of 15 σ fit to the background region. These steps are followed by a standard ramp fitting routine. ExoTiC-JEDI is also able to perform custom bias subtraction, but we find that it does not improve the precision of the data in this case. We extract our Stage 2 products, the 2D wavelength array, and exposure times using the standard jwst pipeline.  \nStage 3 of the ExoTiC-JEDI package performs pixel corrections, additional background and 1/f noise removal, and spectral extraction. Using the data quality  \nflags provided from the jwst pipeline, we replace pixels identified as bad, saturated, low quantum yield, hot, dead, or no gain with the median of the surrounding pixels. To remove additional bad pixels due to cosmic rays or other phenomena, we identify both spatial and timeseries outliers in the data cube with a 20 σ threshold in time and 6 σ threshold spatially, replacing any identified pixels with the median of that pixel in the surrounding 10 integrations or 20 pixels in that row. Any remaining 1/f noise is removed by masking the illuminated region of the detector to calculate the median illuminated pixel value in each column. To extract the 1D stellar spectrum, we fit a Gaussian to each column of the data followed by a fourth-order polynomial to the trace center and widths ( ∼ 0.7 pixels wide). The trace centers and widths are then smoothed with a median filter and used to determine a simple aperture region 5 × the trace FWHM ( ≈ 7 pixels). We use intrapixel extraction to obtain our 1D stellar spectrum. At Stage 3, we also measure the trace position movement on the detector in the x- and y-position for detrending at later stages.  \nWe perform light curve fitting on broadband NRS1 and NRS2 spectra, as well as spectroscopically across the full wavelength range. Using the broadband spectra for each detector and visit, we fit for the planetary system inclination and a/R ⋆ while fixing the period (1.628931 days) and eccentricity (0.0) to literature values presented by Bonfils et al. (2018); as the eccentricity value presented in (Bonfils et al. 2018) is an upper limit we test the fit fixed at both e=0.0 and e=0.22 and find no impact on the resultant transmission spectrum. These parameters, along with the center of transit time, are held constant in the spectroscopic light curve analysis. Stellar limb-darkening coefficients are calculated using the ExoTiC-LD package with a custom model input using a Phoenix stellar model (Husser et al. 2013, Teff=3300 K, logg=5.0, [Fe / H]=0.0) and the nonlinear limb-darkening law. Limb darkening values are then fixed in our light curve analysis. We use a leastsquares optimizer with a batman (Kreidberg 2015) transit model to fit for the transit depth in each bin. We simultaneously fit a series of systematic models to the data and determine the optimal model based on the negative log-likelihood, which incorporates a penalization in complexity based on the AIC (Akaike Information Criterion). We find that the best systematic model, S ( λ ), corrects for a linear trend in time, t , plus the change in x-position, x s , multiplied by the absolute magnitude of the y-positional change, | y s | , such that S ( λ ) = s 0+( s 1 × x s | y s | ) + ( s 2 × t ) , (where s 0 , s 1 , s 2 are coefficient terms).', '3.4. Agreement between Pipelines and Resolutions': 'Figure 1 shows the R ∼ 100 transmission spectra for all three reductions and both visits, demonstrating excellent agreement between the pipelines. We find that our three R ∼ 100 reductions agree better than the native resolution light curve fits, particularly at the red end of the spectrum where the signal-to-noise (SNR) decreases due to a combination of decreased throughput and decreased stellar signal (not shown here). While previous work has found that noise is reduced by performing light curve fits at native pixel resolution prior to spectrally binning the data (Espinoza et al. 2023), our results suggest that low-SNR transit signals can impart biases that can change the shape of the spectrum when fit at native pixel resolution. Lustig-Yaeger & Fu et al. (2023) showed that the precision improvement from fitting at the native pixel resolution can be negated by sufficiently removing 1/f noise at the group level, suggesting that binning light curves prior to fitting is acceptable. Our results here show that binning light curves prior to fitting may be preferable for lowSNR targets (specifically when the spectrophotometric scatter approaches the transit depth). We further discuss the impact of native resolution light curve fitting in Section A.1.', '4. INTERPRETATION': 'To interpret our GJ 1132 b transmission spectra, we first explore the consistency between visits then conduct stellar and planetary atmosphere forward modeling and retrievals on both visits independently.', '4.1. Differences Between Visits': "While the data reduction pipelines show good agreement, the transmission spectra show notable differences between Visit 1 and Visit 2 (see Figure 1). We begin our analysis by investigating the statistical significance of these differences against the null hypothesis of a flat, featureless transmission spectrum. We approach this in two ways, as described below.  \nFirst, we assess how expected or unexpected our measurements would be under the assumption that the spectrum is featureless. We report in Table 1 the reduced chi-squared, χ 2 ν , for the featureless spectrum (flat line model) that best fits each dataset. We then calculate the distribution of expected χ 2 ν , under the assumption that the null hypothesis is true, using 100,000 randomly generated synthetic featureless spectra with the same uncertainties as the observed data. We calculate the probability that χ 2 ν would be at least as extreme as the observed  \nvalue under the assumption that the null hypothesis is correct (i.e., the 'p-value'). Table 1 reports the p-values and corresponding 'sigma' significance with which the null hypothesis is disfavored by the test.  \nSecond, following Moran & Stevenson et al. (2023), we fit each spectrum with a Gaussian model (representing an agnostic spectral feature) and compare it to the featureless model. Both the Gaussian and featureless models are fit using the Dynesty nested sampling code (Speagle 2020), which returns the Bayesian evidence for each fit. From the evidence, we calculate a Bayes factor and convert it into a 'sigma' value (Trotta 2008) representing the significance of the Gaussian feature model over the featureless model. These results are reported in Table 1, where positive values denote evidence favoring the Gaussian model and negative values favor the featureless model. We note that these results can differ from the first statistical test, depending on how well a single Gaussian with varying wavelength center, width, and amplitude can actually fit the spectrum.  \nThe conclusions from our two statistical tests are consistent, despite having slightly different numerical results. In general, Visit 1 contains marginal evidence to reject the null hypothesis and favor a non-flat spectrum, while Visit 2 is statistically consistent with a flat line. For Visit 1, the FIREFLy reduction has a χ 2 ν near unity (indicating it favors the null hypothesis), which disagrees with that of Eureka! and ExoTiC-JEDI -likely owing to the slightly larger uncertainties in the FIREFLy reduction. However, all three reductions agree well for Visit 2 and favor a featureless spectrum.  \nWe also investigated the sensitivity of our results to a potential transit depth offset between the NRS1 and NRS2 detectors for NIRSpec G395H. The second set of rows in Table 1 show the results for the same tests as the top three rows, but now allowing for NRS2 to shift vertically relative to NRS1 to account for a potential systematic offset between the two detectors. These results show that including an offset for NRS2 erodes the statistical significance with which to reject the null hypothesis for Visit 1, while leaving the results for Visit 2 largely unchanged. In Appendix B, Figure 8, we provide a visual example of these null hypothesis tests. These results demonstrate that the two visits are inconsistent unless we allow for a transit depth offset between the two detectors. There is, however, no strong evidence for the necessity of a significant detector offset, since there was no need for a superbias correction step (unlike for GJ 486 b, where this was the primary driver for a detector offset, see Moran & Stevenson et al. 2023). Even when including an offset, Visit 1 is still more consistent with spectral features than Visit 2.  \nTable 1. Is it Flat?  \nNote -χ 2 ν is the reduced chi-squared resulting from the best fitting featureless fit (null hypothesis) to the observed spectrum, 'p-value' 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 (negative values denote preference for the featureless model).", '4.2. Evidence for Variable Star Heterogeneity?': 'A natural physical explanation for different transmission spectra for a planet orbiting an M dwarf is stellar variability (e.g., Rackham et al. 2018). To investigate the possibility of changing starspot coverage, we first perform forward modeling of the out-of-transit extracted, flux-calibrated stellar spectra for both visits.  \nWe compute multi-component stellar forward models using the Allard et al. (2012) PHOENIX models. We seek plausible evidence of evolution or rotation of features onto or off of the observed disk by quantifying the spot and faculae covering fractions of the stellar surface for Visits 1 and 2. We followed a similar procedure to Moran & Stevenson et al. (2023), employing 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 + 100K). We assume that all spots have a common T eff , log( g ), and metallicity (as do the faculae), and constrain each feature to not exceed 45% of the stellar surface. The grid of PHOENIX models used for our analysis covers T eff = 2500-4500 K, log( g ) = 4-5.5 cm s -2 , and [Fe/H] = -0.50, which provides extensive coverage of possible spot and faculae temperatures for GJ 1132. Finally, we assume photospheric values near the literature quoted T eff = 3270 ± 140 K; (Bonfils et al. 2018), log( g ) = 4.88 ± 0.07 cm s -2 (Southworth et al. 2017), and [Fe/H] = -0.12 ± 0.15 (Berta-Thompson et al. 2015).  \nTo compare with the observed baseline spectra, we first convert the native model wavelengths ( ˚ A) and flux densities (erg s -1 cm -2 cm -1 ) to µ m and mJy. We also scaled the models by R 2 ∗ / d 2 using literature values for GJ 1132: R ⋆ =0.21R ⊙ (Bonfils et al. 2018) and d =  \n12.61 pc (Gaia Collaboration et al. 2021). We smoothed and interpolated the models the same resolution as the observations before calculating a χ 2 ν . In our χ 2 ν calculations, we considered 3206 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.1.  \nOur out-of-transit stellar spectra and best-fitting models for both visits are shown in Figure 2. The preferred models for Visit 1 and Visit 2 both have a background photosphere with T eff =3200K, log( g ) = 4.5 cgs, and [Fe/H] = 0, spots with T eff = 2900K, and faculae with T eff = 3500K (top panels, Figure 2). The preferred model for Visit 1 is 33% photosphere, 40% spots, and 27% faculae ( χ 2 ν = 1 . 22). The preferred model for Visit 2 is 35% photosphere, 39% spots and 26% faculae (with a χ 2 ν of 1.16). The largest spectral deviations between visits occur from 2.75-3.3 µ m and 4.3-5.3 µ m for both the observations and the models (bottom panel, Figure 2). However, a 1% difference in spot and faculae coverage is negligible, considering general model uncertainties and error bars on the observations. Therefore, it is unlikely that the differences between visits are caused by the evolution of surface features or rotation onto or off of the visible disk.  \nWe further note that there is no evidence for occulted starspots in either of the visits, which could otherwise explain the different transmission spectra due to the impacts on the morphology of the transit light curves themselves. Figure 6 in Appendix A.2 shows the Eureka! white light curves for both visits, showing the lack of ob-  \nFigure 2. Out-of-transit stellar spectra of GJ 1132 compared to heterogeneous stellar models . Top panel: the extracted GJ 1132 spectra from Visit 1 (blue errors) and Visit 2 (black errors) are compared to best-fitting multi-component PHOENIX models (yellow for Visit 1; pink for Visit 2). The best-fitting Visit 1 model has 40% spot coverage and 27% faculae coverage, compared to 39% spot coverage and 26% faculae coverage for Visit 2. Middle panels: zoom-ins on the gray highlighted regions in the top panel. Bottom panel: deviations between visits for the models (yellow) and observations (black). The largest deviations for both models and observations occur near 3 µ m and long-ward of 4.3 µ m.  \n<!-- image -->  \nvious spot occultations in either visit (occulted starspots would result in a brief decrease in transit depth during the transit event). However, because occulted starspots typically have a minimal impact on the visible shape of a light curve at these wavelengths, we also consider how the fitted orbital parameters may change due to such features. FIREFLy and ExoTiC-JEDI fit the white light curves from the two visits independently, and while the independent fit values are not reported here, they are consistent within uncertainties.', '4.3. An Atmosphere Around GJ 1132 b?': "We next assess possible atmospheric explanations for GJ 1132 b's transmission spectrum. As in our previous studies (Lustig-Yaeger & Fu et al. 2023; Moran & Stevenson et al. 2023), we compare each reduction to a set of simple forward models to compare possible atmospheric compositions and a no-atmosphere sce-  \nnario. We first generate forward model atmospheres using either thermochemical equilibrium CHIMERA (Line & Yung 2013; Line et al. 2014) models or simple one-or-two gas isothermal atmosphere PICASO (Batalha et al. 2019) models. We then compute model transmission spectra using PICASO 's radiative transfer module. We bin the resulting model spectra to the resolution of each reduction and compute a χ 2 ν (with 60 degrees of freedom (dof) for the Eureka! and FIREFLy reductions and 59 for the ExoTiC-JEDI reduction) to assess goodness of fit. We summarize our results in Figure 3.  \nFirst, we run a set of thermochemical equilibrium models with CHIMERA . We include opacities from H, collision-induced absorption (CIA), H 2 , He, H 2 O, CH 4 , CO, CO 2 , NH 3 , N 2 , HCN, and H 2 S and use the parameterized temperature-pressure profile of Guillot (2010) with an equilibrium temperature of 530 K. We run these forward models at 100 × to 1000 × solar metallicities,  \nFigure 3. Transmission spectra of GJ 1132 b compared to atmospheric forward models. The Eureka! reduction (black circles) is shown for both visits at an R ∼ 100, with Visit 1 on the left and Visit 2 on the right. Also shown are a series of end-member atmospheric forward models generated with PICASO to illustrate goodness of fit ( χ 2 ν ) of various scenarios for each visit, shown after each model for Visit 1 and 2, respectively. Dashed lines show poorer overall fits including a 1000 × solar atmosphere (dashed purple), a pure 1 bar CH 4 atmosphere (dashed navy), and a pure 1 bar CO 2 atmosphere (dashed orange). Better fits include pure H 2 O atmospheres (blue) or water atmospheres with 10% N 2 O (yellow) or 10% CH 4 (pink), or a flat line indicative of no atmosphere or a high altitude opaque aerosol layer (gray dashed). In Visit 1, water-rich atmospheres with other species are preferred, but the features driving this fit disappear in Visit 2 so that an atmosphere-free model is the best fit.  \n<!-- image -->  \nfinding that the scale height of the atmosphere is so large in all cases that we obtain poor fits for all reductions in Visit 1 ( χ 2 ν ≳ 1.43). Therefore, we rule out clear hydrogen-dominated atmospheres in thermochemical equilibrium at moderate confidence ( ≳ 2.5 σ ). For Visit 2, these confidences are decreased for the 1000 × solar metallicity case (down to 1.3 σ for the FIREFLy reduction), but a hydrogen-dominated atmosphere is never the statistically preferred scenario compared to our forward models. A long-lived hydrogen-dominated atmosphere would be unexpected given the planet's density and radius (Luger et al. 2015; Rogers 2015; Estrela et al. 2020; Rogers et al. 2021), so our forward model limits here fit expectations better than some previous results for GJ 1132 b (Southworth et al. 2017; Swain et al. 2021).  \nFor our simpler, non-self-consistent PICASO models, we examine whether 1bar, isothermal atmospheres of pure CH 4 , pure CO 2 , or pure H 2 O agree with the data from each reduction. We find that CH 4 -dominated atmospheres (dark blue dashed line in Figure 3) are the most strongly ruled out for both visits, to at least 4.2 σ across all reductions. Although in Visit 1 there is a rise in transit depth at the strong CH 4 absorption feature centered at 3.3 µ m, the lack of strong CH 4 absorption  \nat the wavelengths probed by NRS2 result in a poor fit. Similarly, CO 2 -dominated atmospheres poorly fit Visit 1. However, H 2 O-rich atmospheres (thick solid lines in Figure 3) provide better fits to Visit 1 due to the broad H 2 O absorption slope at the bluest wavelengths probed by NIRSpec G395H. In Visit 1, an uptick in transit depth at ∼ 4.5 µ mis also noticeable. Multiple molecules, including O 3 , CS 2 , and N 2 O, have an absorption band around this wavelength (e.g., Schwieterman et al. 2022), but of these N 2 O has the best-matching feature center and width. Therefore, in addition to our pure H 2 O atmosphere, we also generate atmospheric models with H 2 O as the background gas and either 10% N 2 O or 10% CH 4 . These result in visually better fits at the two increases in transit depth at 3.3 µ mand 4.5 µ m. A flat-line fit is slightly disfavored for Visit 1 (rejected between 13 σ ) - consistent with our previous statistical tests in Section 4.1 - in favor of H 2 O-rich atmospheres.  \nIn summary, our forward models prefer a H 2 Odominated atmosphere for Visit 1. While adding 10% CH 4 or N 2 O visually explain the observed features, they do not improve χ 2 ν (but we only consider a single mixing ratio for both species). We explore the full range of possible mixing ratios consistent with Visit 1 with atmospheric retrievals in Section 4.4.  \nFor Visit 2, we find that a flat line - indicative of either no atmosphere or a high altitude aerosol layer produces the lowest χ 2 ν (consistent with Section 4.1). At GJ 1132 b's ∼ 500 K equilibrium temperature, while condensate clouds are unlikely (given the lack of condensible species) photochemical hazes could form in a variety of atmospheres (e.g., Horst et al. 2018; He et al. 2018; Gao et al. 2020). Moreover, we cannot rule out the H 2 O-dominated atmosphere preferred by Visit 1 from the Visit 2 data, given its low χ 2 ν ( < 1).", '4.4. Retrieval Analysis: An Atmosphere or Starspots?': "Our analysis thus far has yielded two key results: (i) GJ 1132's out-of-transit stellar spectrum strongly favors a heterogeneous star with constant spot and faculae properties between the two Visits, and (ii) GJ 1132 b's Visit 1 spectrum can be explained by an atmosphere, but Visit 2 is statistically flat. Here we attempt to reconcile these insights via retrieval modeling of GJ 1132 b's transmission spectrum considering both atmospheric and unocculted starspot scenarios. Our retrieval results are summarized in Figure 4.", '4.4.1. Atmosphere Scenario': "We first explore the range of atmospheres consistent with GJ 1132 b's transmission spectrum via separate retrievals of our two visits using the open source POSEIDON code (MacDonald & Madhusudhan 2017; MacDonald 2023). We considered 11 potential gases to span a wide parameter space of plausible atmospheric compositions: N 2 , H 2 , H 2 O, CO 2 , CH 4 , N 2 O, NO 2 , HCN, NH 3 , SO 2 , and PH 3 . The opacities used for the retrieval forward model are described in MacDonald & Lewis (2022). The mixing ratios of these gases can range from 10 -12 to 1, using centered log-ratio priors as in Lustig-Yaeger & Fu et al. (2023) and Moran & Stevenson et al. (2023). The other free parameters are (priors in brackets) the atmospheric temperature ( U (400 K, 900 K)), the atmosphere radius at the 10 bar reference pressure ( U (0.85 R p , obs , 1.15 R p , obs )), the haze power-law exponent ( U (-20, 2)) and log-Rayleigh enhancement factor ( U (-4, 8)) - defined as in MacDonald & Madhusudhan (2017) - and the log-pressure of an opaque cloud/surface ( U (-7, 2), in bar). We calculate transmission spectra via opacity sampling at a resolving power of R = 20,000 from 0.65.2 µ m, before convolving the model with the instrument point spread function and binning to the resolution of the observations. We sample this 15-parameter space using the PyMultiNest (Feroz et al. 2009; Buchner et al. 2014) package with 2,000 live points.  \nWe perform retrievals on each GJ 1132 b visit separately to consider how their different morphology (discussed in Section 4.1) affects atmospheric inferences.  \nOur retrievals here focus on the Eureka! reduction since we found consistent retrieval results with ExoTiC-JEDI and FIREFLy . We do not consider a free detector offset between the NRS1 and NRS2 spectra in these retrievals, given the lack of evidence for a superbias correction during the data reduction.  \nOur Visit 1 retrieval (left panels of Figure 4) favors a H 2 O-dominated atmosphere (Bayes factor = 6.5 / 2.5 σ ) with trace amounts of CH 4 (Bayes factor = 27 / 3.1 σ ). Figure 4 demonstrates that the evidence for H 2 O is driven by the slope seen in the NRS1 data, while a feature near 3.3 µ m is attributed to CH 4 . A weak feature near 4.5 µ m is best fit by N 2 O absorption, but our Bayesian model comparison indicates insufficient evidence for N 2 O (Bayes factor = 1.2 / 1.4 σ ). The volume mixing ratio posterior distributions in Figure 4 show a H 2 O abundance consistent with 100% (log X H 2 O = -0 . 01 +0 . 01 -0 . 43 ), a CH 4 abundance of ∼ 1% (log X CH 4 = -2 . 37 +1 . 13 -0 . 69 ), and an N 2 O abundance of ∼ 100 ppm (log X N 2 O = -4 . 64 +1 . 54 -4 . 24 ). We note that our Visit 1 retrieval finds no evidence of HCN, which was previously suggested from the HST WFC3 analysis of GJ 1132 b by Swain et al. (2021). However, our G395H data does not rule out a scattering slope below 2 µ m (as shown by the wide 2 σ confidence region in Figure 4) similar to that inferred by Swain et al. (2021) (but cf. Libby-Roberts et al. 2022 and Mugnai et al. 2021).  \nIn contrast, our Visit 2 retrieval (right panels of Figure 4) is consistent with a flat line with no constraints on the atmospheric composition. Such a flat spectrum can be explained by many degenerate atmospheric properties, including a high mean molecular weight, low temperature, low surface pressure, and/or a high-altitude aerosol layer. Given that our data are sufficiently precise to differentiate between several high mean molecular weight atmospheres with high surface pressure / deep clouds (as shown by our inference of a H 2 O-dominated atmosphere from Visit 1), our favored explanation for the featureless Visit 2 spectrum is a high-altitude aerosol layer. However, a wide range of cloud top pressures are permitted (log P cloud < -1 . 3 to 1 σ ; see Figure 4) after marginalization over all the possible combinations of atmospheric temperature and background gases with higher mean molecular weight than H 2 O.  \nAssuming GJ 1132 b's transmission spectra are explained by a planetary atmosphere, our retrievals thus suggest that the cloud opacity would need to significantly increase between Visit 1 and 2 to explain our different spectra. As further discussed in Section 4.5, it is highly improbable that an atmosphere can change from a relatively clear state in one visit to host a global high-altitude cloud in the next visit. Furthermore,  \nFigure 4. Atmospheric and starspot retrieval results for GJ 1132 b. Top panels: comparison between the retrieved transmission spectra for the Eureka! reduction of Visit 1 (left) and Visit 2 (right) adopting two distinct retrieval models: (i) a planetary atmosphere with no unocculted starspots (green contours), and (ii) no atmosphere with unocculted starspots (orange contours). The median retrieved spectrum (solid lines) and 1 σ and 2 σ model confidence intervals (dark and light contours) for each scenario are overlaid. Labels indicate the locations of H 2 O, CH 4 , and N 2 O absorption bands. Middle panels: posterior histograms for the atmosphere scenario, highlighting the volume mixing ratios of the three molecules tentatively inferred from Visit 1: H 2 O, CH 4 , and N 2 O. Bottom panels: posterior histograms for the unocculted starspot scenario, defined by the spot coverage fraction, spot temperature, and the background stellar photosphere temperature. Visit 1 requires either a water-rich atmosphere with trace CH 4 (and possibly N 2 O) or unocculted starspots, but the scenarios cannot be differentiated without observations shortwards of 3 µ m. Visit 2 is, however, consistent with no atmosphere and no unocculted starspots.  \n<!-- image -->  \nGJ 1132 b's equilibrium temperature places it in a parameter space without obvious condensable material for clouds to form (e.g., Gao et al. 2021), which would suggest a photochemical haze as the cause of this aerosol layer. Similarly, a transition from a relatively clear atmosphere to one with such high haze opacity would be highly improbable given no change in radiative forcing. However, given the wide uncertainty on our retrieved cloud-top pressure, we note that both Visits 1 and 2 are consistent with an intermediate cloud pressure at ∼ 10 mbar to 1 σ . We next turn to consider an alternative explanation that does not require an atmosphere.", '4.4.2. Unocculted Starspot Scenario': 'We next investigate the alternative explanation of unocculted starspots shaping the observed transmission spectra, adopting the same retrieval configuration as Moran & Stevenson et al. (2023). This retrieval model assumes no planetary atmosphere, with unocculted stel-  \nlar heterogeneities producing any wavelength-dependent features in the transmission spectrum (see Rackham et al. 2023, for a review of the impact of starspots on transmission spectra). Four parameters define this model (priors in brackets): the heterogeneity temperature, T het ( U (2300 K, 1.2 T ∗ , eff )), the heterogeneity covering fraction, f het ( U (0.0, 0.6)), the photospheric temperature, T phot ( N ( T ∗ , eff , σ T ∗ , eff )), and the planetary radius, R p ( U (0.9 R p , obs , 1.1 R p , obs )). The priors are specified in terms of literature properties of the host star: T ∗ , eff = 3270K and σ T ∗ , eff = 140K (Bonfils et al. 2018). We calculate the impact of the transit light source effect (Rackham et al. 2018) by interpolating PHOENIX models (Husser et al. 2013) using the PyMSG package (Townsend & Lopez 2023). We verified that a more complex parameterization of stellar contamination, including both faculae and spots, does not improve the fit or alter the retrieved spot properties, so we focus here on results assuming a single-heterogeneity population.  \nAs in the previous section, we perform independent retrievals for the two visits using the Eureka! reduction and without a free offset between the detectors.  \nOur Visit 1 retrieval is well-explained by unocculted starspots covering ∼ 20% of the stellar surface with a spot temperature ∼ 400 K cooler than the photosphere. These starspot properties are consistent with the stellar spectrum fits described in Section 4.2. As shown in Figure 4, the starspot scenario explains the Visit 1 data with a spectral slope shortwards of 4 µ m. The posterior distributions in Figure 4 demonstrate that a wide range of spot coverage fractions are consistent with Visit 1 ( f het = 0 . 23 +0 . 21 -0 . 10 ) - due to the f het -T het degeneracy (e.g. Rathcke et al. 2021, their Figure 10).  \nHowever, once again our Visit 2 retrieval is consistent with a flat, featureless spectrum. Under the starspot model, this requires either a low spot coverage fraction or a spot temperature similar to the stellar photosphere. However, the retrieved starspot fraction for Visit 2 is formally consistent with Visit 1 within 1 σ , which agrees with the out-of-transit stellar modeling in Section 4.2 that found a negligible difference between visits.', '4.4.3. An Atmosphere vs. Starspots': "At first glance, the available statistical evidence equally supports the atmosphere scenario or the starspot scenario for Visit 1. The Bayesian evidences (ln Z = 490 . 1 and 490 . 3 for the atmosphere vs. starspot models, respectively) and minimum reduced χ 2 ( χ 2 ν = 1 . 35 and 1 . 36 for the atmosphere vs. starspot models, respectively) are indistinguishable. However, this is largely a manifestation of the higher dimensionality of our reference atmosphere model compared to the starspot model (15 parameters vs. 4 parameters for the starspot model). Comparing the χ 2 directly, we find the maximum likelihood atmospheric scenario model ( χ 2 = 61 with 45 dof) achieves a better quality fit than the starspot scenario ( χ 2 = 76 with 56 dof). This better fit arises from the spectral slope caused by H 2 O absorption in combination with the CH 4 feature near 3.3 µ m providing a better fit than the slope produced by unocculted starspots.  \nWe, therefore, ran an additional, simplified, atmosphere scenario retrieval focusing only on those properties indicated by the data. Our motivation is to account for the 'Occam penalty' disfavoring models with redundant free parameters (i.e., the non-detected molecules in the atmosphere scenario). We thus consider a H 2 Odominated atmosphere with trace CH 4 and N 2 O, constituting a 4-parameter model defined by the planetary reference radius, temperature, and the CH 4 and N 2 O abundances (the H 2 O abundance determined via abundances summation). This simplified model obtains an  \nexcellent fit to Visit 1 ( χ 2 = 63 with 56 dof; χ 2 ν = 1 . 13; ln Z = 494 . 2), providing tentative evidence favoring an atmosphere over starspots for Visit 1. The lack of observed starspot variability in the out-of-transit stellar spectrum (Figure 2), in small tension with the low spot fraction required to render Visit 2 flat (Figure 4), may also support the atmosphere interpretation.  \nHowever, the inferred atmosphere from Visit 1 - a H 2 O-dominated atmosphere with trace CH 4 and N 2 O - would be unstable against a runaway greenhouse, as all three molecules are highly effective greenhouse gases and highly susceptible to photolysis (e.g., Rugheimer et al. 2015). Thus, such an atmosphere would rapidly be lost to space (e.g., Goldblatt et al. 2013), requiring ongoing outgassing, a very high ( > 5 wt%) initial H 2 O inventory at formation, and a present-day magma ocean (Schaefer et al. 2016) to replenish the H 2 O to the levels suggested by our Visit 1 atmospheric retrieval. Moreover, an H 2 O-dominated atmosphere would be expected to also include ample oxidized carbon species, such as CO 2 or CO, given outgassed abundances suggested from previous studies (Sossi et al. 2023; Tian & Heng 2023). While we do not see evidence for CO 2 , our NIRSpec G395H data cannot constrain the presence or absence of CO. The co-presence of ∼ percent levels of CH 4 with CO in an H 2 O-dominated atmosphere would require finely tuned carbon abundances and oxygen fugacities of the planetary interior (Tian & Heng 2023).  \nGiven the disagreement between our two visits, we suggest it is premature to claim a clear preference for an atmosphere or starspots. Nevertheless, our JWST observations do rule out the H 2 -dominated atmosphere with HCN previously suggested by Swain et al. (2021) (see Appendix C). One path to a resolution is the significant predicted deviations between the atmosphere and starspot scenarios at shorter wavelengths (see Figure 4), which could be probed with future observations.", '4.5. Potential For Planetary Variability': "The impact of atmospheric variability on transmission spectra of tidally locked terrestrial planets has been explored by numerous teams (e.g. May et al. 2021; Song & Yang 2021; Cohen et al. 2022; Rotman et al. 2023), with the general conclusion being that such variability is below the detectability limit of JWST instruments. This is because, while the cloud cover of any one region of the terminator can be highly variable, we are observing limb averaged spectra that simultaneously probe cloudy and cloud-free regions. With the difference between GJ 1132 b's two spectra of order ∼ 100 ppm at the 4.5 µ m 'feature' (see Figure 1), any planetary variability would be required to be at least an order of magni-  \ne larger than predicted by current GCMs of temperate, tidally locked terrestrial planets in the above works. While GJ 1132 b is hotter than the models in the above works, it remains physically unlikely that such variability is the cause of the observed spectral differences. This is because at the equilibrium temperature of GJ 1132 b the planetary dayside is more likely to be cloud free. Further, such large scale variability in transmission has yet to be detected on even the most optimal targets hot Jupiters (e.g. Kilpatrick et al. 2020).", '5. CONCLUSIONS': "We presented two transit observations of the superEarth GJ 1132 b with JWST NIRSpec G395H that yield distinctly different transmission spectra. Barring other possibilities, the differences between the visits may be explained by random noise fluctuations that took the unfortunate shape of spectral features for one visit. Without a third observation of GJ 1132 b at these wavelengths, it is impossible to determine which of the two visits most accurately reflects the true nature of the planet. Our conflicting observations demonstrate the potential risk of claiming atmospheric detections for rocky exoplanets based on a single JWST observation.  \nShould Visit 1 represent the 'truth' for GJ 1132 b our retrievals exhibit a slight preference for a H 2 Odominated atmosphere with trace CH 4 and N 2 O ( χ 2 ν = 1 . 13) compared to contamination from unocculted starspots ( χ 2 ν = 1 . 36). These two scenarios would produce significantly different transmission spectra at shorter wavelengths (see Figure 4). While HST WFC3 and ground-based optical data exist for GJ 1132 b, their lack of wavelength overlap with our NIRSpec G395H observations results in too much freedom when accounting for a possible inter-instrument transit depth offset, resulting in either scenario being allowed (see Appendix C). However, our Visit 1 observations do rule out a H 2 -dominated atmosphere containing HCN, as previously suggested by the HST WFC3 analysis of Swain et al. (2021). See Appendix C for further comparison of our new JWST observations to the existing HST WFC3 and ground-based data.  \nInstrument systematics could provide an alternative explanation for our divergent transmission spectra. Allowing for a possible offset between the NRS1 and NRS2 detectors would eliminate the statistical significance of the bluewards slope in the Visit 1 data (see Section 4.1 and Appendix B), which drives our inference of an atmosphere or unocculted starspots. Larger detector offsets were seen between NRS1 and NRS2 in Moran & Stevenson et al. (2023), but were manually corrected based on differences between the visits and not explored  \nin their equivalent Gaussian feature tests. However, we note that we do not see evidence during our data reduction process for the need for a similar superbias correction, and find equivalent spectra regardless of applying such a correction. Future observations of GJ 1132 b with NIRSpec G395M would cover the same wavelengths without the detector gap - to determine the true shape of GJ 1132 b's transmission spectrum at these wavelengths. With the growing evidence for detector offsets within NIRSpec G395H data, it may be preferable to observe planets with small predicted spectral features using NIRSpec G395M instead, especially if the star is dim enough to allow it, to ensure that detector offsets are not mistaken for atmospheric features.  \nAdditional future observations at wavelengths below 3 µ m, but still overlapping with the NIRSpec G395H wavelength range, are crucial to breaking the degeneracy between the atmospheric and unocculted starspot explanations. This atmosphere-starspot degeneracy has now been seen in multiple NIRSpec G395H observations of small planets orbiting M dwarfs (see also Moran & Stevenson et al., 2023). For example, NIRISS SOSS observations combined with our existing NIRSpec G395H data, or new NIRSpec G395M data (avoiding any potential detector offsets), will be crucial to determine if (1) GJ 1132 b's spectral features are real and repeatable, and (2) if these features are best described by a planetary atmosphere or stellar contamination.  \nWhile the quest continues for an unambiguous atmospheric detection on a rocky planet, our results for GJ 1132 b provide an important reminder of the necessity of repeat observations to confirm the reliability of potential detections.", 'ACKNOWLEDGMENTS': "This work is based in part on observations made with the NASA/ESA/CSA JWST . The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. 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. Support for JWST program #1981 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127. We thank the anonymous referee for their constructive and timely feedback, alongside fruitful discussions  \nthat consequently improved our study. R.J.M. is supported by NASA through the NASA Hubble Fellowship grant HST-HF2-51513.001, also awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 5-26555. Work by S.P. is supported by NASA under award number 80GSFC21M0002. J.K. acknowledges financial support from Imperial College London through an Imperial College Research Fellowship grant. This material is based upon work performed as part of NASA's CHAMPs team, supported by the National Aeronautics and Space Administration (NASA) under Grant No. 80NSSC21K0905 issued through the Interdisciplinary Consortia for Astrobiology Research (ICAR) program. The authors thank Raissa Estrela and the Excalibur team at JPL for useful discussions on the HST model presented in Swain et al. (2021).  \nFacilities:  \nJWST (NIRSpec)  \nThe JWST data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute. The spe-  \ncific observations analyzed can be accessed via DOI: 10.17909/0njr-8110.", '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-JEDI (Alderson et al. 2022), ExoTiC-LD (Grant & Wakeford 2022), FIREFLy (Rustamkulov et al. 2022) IPython (P'erez & Granger 2007), jwst (Bushouse et al. 2022), lacosmic (van Dokkum 2001), lmfit (Newville et al. 2014), Matplotlib (Hunter 2007), NumPy (van der Walt et al. 2011; Harris et al. 2020), PHOENIX (Allard et al. 2012), PICASO (Batalha et al. 2019), POSEIDON (MacDonald & Madhusudhan 2017; MacDonald 2023),  \nPyMultiNest (Feroz et al. 2009; Buchner et al. 2014), SciPy (Virtanen et al. 2020)", 'REFERENCES': "Birkmann, S. M., Ferruit, P., Giardino, G., et al. 2022,  \nA&A, 661, A83, doi: 10.1051/0004-6361/202142592  \nBonfils, X., Almenara, J. M., Cloutier, R., et al. 2018,  \nA&A, 618, A142, doi: 10.1051/0004-6361/201731884  \nBourque, M., Espinoza, N., Filippazzo, J., et al. 2021, The Exoplanet Characterization Toolkit (ExoCTK), 1.0.0, Zenodo, doi: 10.5281/zenodo.4556063  \nBuchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125, doi: 10.1051/0004-6361/201322971  \nBushouse, H., Eisenhamer, J., Dencheva, N., et al. 2022, JWST Calibration Pipeline, 1.8.2, Zenodo,  \ndoi: 10.5281/zenodo.7325378  \nCohen, M., Bollasina, M. A., Palmer, P. I., et al. 2022, ApJ, 930, 152, doi: 10.3847/1538-4357/ac625d  \nDiamond-Lowe, H., Berta-Thompson, Z., Charbonneau, D., & Kempton, E. M. R. 2018, AJ, 156, 42,  \ndoi: 10.3847/1538-3881/aac6dd  \n- do Amaral, L. N. R., Barnes, R., Segura, A., & Luger, R. 2022, ApJ, 928, 12, doi: 10.3847/1538-4357/ac53af  \nDong, C., Jin, M., & Lingam, M. 2020, ApJL, 896, L24, doi: 10.3847/2041-8213/ab982f  \nEspinoza, N., ' Ubeda, L., Birkmann, S. M., et al. 2023, PASP, 135, 018002, doi: 10.1088/1538-3873/aca3d3", 'A. DATA ANALYSIS': 'A.1. Light Curve Fitting: To Spectrally Bin Before or After?  \nLight curves from Eureka! and FIREFLy were extracted at both the R ∼ 100 resolution and at the native pixel level resolution. The native pixel level resolution light curves were fit in the same manner as described in Sections 3.1 and 3.2 and then binned to the same wavelength grid as the R ∼ 100 spectra to compare the effect of binning before and after light curve fitting. Figure 5 shows the differences between Eureka! and FIREFLy transmission spectra with these two methods. We find that the transmission spectra diverge significantly on the red side of the NRS2 detector when fitting at the native pixel resolution and then spectrally binning, likely due to the lower SNR at longer wavelengths. This suggests that fitting at the native resolution is not always appropriate.  \nNotably, it is the Eureka! reduction that changes the most when binning before or after light curve fitting. We conducted a further comparison using a least squared fitting routing in Eureka! (instead of emcee ) and find better agreement with FIREFLy in the native pixel light curve fits. However, an analysis of the emcee chains for both resolutions shows no evidence that either fit has not converged and no multi-modal distributions are present in the posteriors. Nevertheless, because the resulting spectra were consistent when binning prior to light curve fitting regardless of the tool used to perform the light curve fitting, we chose to adopt the pre-binned R ∼ 100 spectra in our transmission spectra analysis.  \nChoices in limb darkening treatment (fixing vs. fitting) have been shown to impact JWST transmission spectra (e.g. Carter & May et al., in prep.) for larger planets. Here we apply their recommendations, specifically to fix limb darkening rather than fit for it, and to pre-bin the data to a lower resolution. This minimizes additional uncertainty and increases accuracy in the final transmission spectrum. We do not see significant differences between the prebinned FIREFLy and Eureka! spectra even with the different treatment of limb darkening between the two reductions ( FIREFLy uses the white light parameters for all spectroscopic bins, Eureka! fixes to wavelength dependent limb darkening based on stellar models). It is therefore unlikely that our limb darkening choices are causing the differences between the pre- and post-binned transmission spectra in this work.  \nFigure 5. Comparison of Eureka! and FIREFLy transmission spectra from binning before and after light curve fitting. In yellow are the Eureka! reductions and in pink are the FIREFLy reductions. For both, the shaded circles are the reductions that are first binned to R ∼ 100 and then fit (pre-binned) and the darker triangles are the reductions that are fit at the native resolution and then binned to R ∼ 100 (post-binned). In the bottom panels we show the differences between the Eureka! and FIREFLy reductions for both cases in units of σ , with the solid thicker line the pre-binned case and the thinner dashed line the post-binned case. The shaded regions denote ± 1 σ and ± 2 σ to guide the eyes. The Eureka! and FIREFLy pre-binned spectra are clearly more consistent with each other (always within a ± 1 σ difference) than the post-binned case, suggesting that fitting at the native resolution is not always appropriate.  \n<!-- image -->', 'A.2. 2D Light Curve Comparisons': 'Figure 6 shows the R ∼ 100 2D light curves for all three reductions for Visit 1. We note that while each analysis has different levels of correlated noise in the data panel (left), the noise is generally well modeled as seen by the models and residuals panels (middle and right). In the bottom panel we show the Eureka! white light curves from Visits 1 and 2 to highlight the precision of the data. The residuals (lower middle) are well described by Gaussians (lower right).  \n<!-- image -->  \nFigure 6. 2D light curves for all reductions from Visit 1 and Eureka! white light curves for both visits . In the top panel we show the 2D light curves for all three reductions, including the data (left), models (middle) and residuals (right) at the R ∼ 100 resolution. FIREFLy trims extra data from the start of the NRS1 data due to a temporal ramp. In the bottom panel we show the Eureka! white light curves and models for NRS1 and NRS2 (left), the residuals (middle), and a histogram of the residuals for both NRS1 and NRS2 (right) for both visits 1 and 2. There is no evidence of occulted starspots in any of the white light curves.  \n<!-- image -->  \nTable 2. Best fit orbital parameters from white light curve fitting.  \nNote -We adopt the Eureka! results as our system parameters. The Eureka! best-fit values are derived from a mean of independent joint-visit fits to the two detectors. Because Eureka! performs a joint fit of the two visits, only one T 0 value is reported. The FIREFLy and ExoTiC-JEDI values are derived from a weighted mean of four independent white light curve fits, one per detector for each visit.', 'A.3. Best-Fit Orbital Parameters': 'In Table 2 we include the best-fit orbital parameters from Eureka! , FIREFLy , and ExoTiC-JEDI . Eureka! performs a joint fit across the two visits independently for the two detectors (a total of two white light fits) while FIREFLy and ExoTiC-JEDI fit all four white light curves independently. The values provided in Table 2 are weighted means of the respective number of white light curve fits from each pipeline.  \nFor spectroscopic light curve fitting, Eureka! and FIREFLy use the values in Table 2 for both visits and detectors while ExoTiC-JEDI uses the best-fit white light parameters from the respective white light fit for each visit and detector.  \nWe note that because the scaled semi-major axis (a/R s ) and orbital inclination (i) are by definition degenerate with one another, all three pipelines converge to a different part of parameter space. Regardless of these differences in a/R s and inclination, no transit depth offset is required to align the spectra from the three reductions (see Figure 1).', 'A.4. Correlated Noise': "All three data reduction pipelines consider correlated noise by inflating error bars, if needed, at the white and spectroscopic light curve stage based on the scatter of the data. The Allan variance plot for Eureka! 's white and spectroscopic light curves (Figure 7) demonstrates that the correlated noise is small enough to not impact our results. For Eureka! , correlated noise in the white light curve appears when binning by timescales greater than 0.5-5 minutes, and this is believed to be due to some inherent thermal cycling from the telescope (Rigby et al. 2022). Because this timescale is much shorter than the transit duration, it is unlikely to affect the white light curve transit depth. Furthermore, correlated noise is negligible in the spectroscopic light curves.  \n<!-- image -->  \nFigure 7. Allan variance plots for the Eureka! reduction. Left panel: variation of the NRS1 and NRS2 detector white light curve precision (root mean square; RMS) with longer bin timescales (more integrations). As seen in previous works (e.g., Moran & Stevenson et al. 2023), correlated noise emerges around timescales of 0.5-5 minutes (5-50 integrations) and is likely due to thermal cycling (Rigby et al. 2022). Right panel: same, for spectroscopic light curves. The individual spectroscopic light curves show no evidence of correlated noise.  \n<!-- image -->", 'B. IS IT FLAT? GAUSSIAN FEATURES AND DETECTOR OFFSETS': "Figure 8 shows an example of the results presented in Section 4.1 and Table 1 that are used to determine the statistical significance of any candidate features in the transmission spectrum. Specifically, Figure 8 highlights the Eureka! reduction for our GJ 1132 b Visit 1 observations. The top row shows the nominal model (a single Gaussian feature), while the bottom row shows the model that allows for an offset between the NRS1 and NRS2 detectors (parameterized by a step function in the model). For both top and bottom, the left panels show the data compared to the best-fitting model, while the right panels show the distribution of expected reduced chi-squared values given the uncertainty on the data points.  \nFor the Eureka! Visit 1 reduction shown, this demonstrates how allowing for the offset between the NRS1 and NRS2 detectors removes the significance of the Gaussian feature detection, and also changes the location of the Gaussian feature detected (from the blue wards slope to the feature near 4.5 µ m). A complete figure set showing analogous results for both visits plus the combined spectrum for each of the three independent reductions is available in the online journal.  \nFig. Set 1. Gaussian Feature Detection TestsFigure 8. Spectral feature detection tests for the Eureka! Visit 1 . The top two panels show baseline results directly fitting the spectrum as shown, while the bottom two panels show results allowing a model offset for the NRS2 detector data ( > 3 . 75 µ m). The left panels show the spectrum compared to the median best-fitting Gaussian feature model with 1 σ and 3 σ credibility envelopes displayed. Positive (negative) 'sigma' values in the legend indicate the confidence with which the Gaussian feature model is favored (disfavored) compared to a featureless model (the null hypothesis). The right panels show the distribution of expected χ 2 ν values assuming the null hypothesis. The observed χ 2 ν is shown as the vertical dashed line along with the percent of the distribution that would be expected to have at least as large of a result. The complete figure set (9 images, for two visits plus the combined spectrum for three independent reductions) is available in the online journal.  \n<!-- image -->", 'C. CONSISTENCY WITH PREVIOUS DATA': 'Here we consider the consistency of our JWST NIRSpec G395H transmission spectra of GJ 1132 b with previously reported near-infrared and optical observations. We first demonstrate that our JWST data rule out the H 2 -dominated atmosphere with HCN and aerosols proposed by Swain et al. (2021). Second, we show that adding the ground-based observations from Diamond-Lowe et al. (2018) is insufficient to differentiate the atmosphere and starspot scenarios.', 'C.1. HST WFC3': "Figure 9 compares our best-fitting retrieved spectrum from the atmosphere scenario (pink line, see Section 4.4.1) to both our JWST data (black circles) and the HST WFC3 reductions presented in Swain et al. (2021) (gray stars) and Libby-Roberts et al. (2022) (black triangles). Our best-fit model for the JWST data is inconsistent with the spectral slope and HCN features present by Swain et al. (2021)'s HST WFC3 reduction. Our best-fitting model is, however, consistent within error (and a modest offset, not applied here) with the reduction presented by Libby-Roberts et al. (2022). We also show an extrapolated variant (purple dashed line) of the Excalibur 3 spectral fit from Swain et al. (2021) (purple solid line) extended from 1.7 to 5 µ mbased on their stated model parameters. This model is inconsistent with the atmospheric retrieval of our observed JWST data, which, as discussed in Section 4.4.1, find no evidence for HCN.  \nFigure 9. Comparison between a proposed H 2 -dominated spectrum of GJ 1132 b and our best-fitting model . The best-fitting H 2 O-dominated atmosphere spectrum from the POSEIDON retrieval (pink curve), which we fit to only the JWST NIRSpec G395H data (black circles), is more consistent with the HST WFC3 reduction from Libby-Roberts et al. (2022) (black triangles) than the reduction from Swain et al. (2021) (gray stars). The best-fit model from Swain et al. (2021) (purple solid curve), extrapolated out to 5 µ m with PICASO (purple dashed curve), is inconsistent with our JWST NIRSpec data.  \n<!-- image -->", 'C.2. Ground-Based Data': 'The ground-based observations of GJ 1132 b presented by Diamond-Lowe et al. (2018) cover important wavelengths (0.7-1.0 µ m) that can potentially differentiate between the atmosphere and starspot scenarios from our JWST NIRSpec G395H Visit 1 retrievals (see Figure 4 and Section 4.4.3). However, due to the different telescopes, instruments, data reduction techniques, and the large wavelength gap between the two datasets, one must consider a free offset in retrievals when including both the ground-based and JWST data.  \nFigure 10 shows our retrieved transmission spectra and atmospheric/stellar properties from two POSEIDON retrievals including the ground-based data from Diamond-Lowe et al. (2018). We find that the retrieved properties are broadly consistent with our retrievals only including our JWST NIRSpec G395H data (see Figure 4). The only differences are a higher preferred CH 4 abundance for the atmosphere scenario and a tighter constraint on the spot coverage fraction for the starspot scenario. However, the inclusion of an offset results allows the ground-based data to move to accommodate either model. Consequently, the ground-based from Diamond-Lowe et al. (2018) does not allow us to differentiate between the starspot and atmosphere scenarios.  \nFigure 10. Retrieval results including ground-based observations of GJ 1132 b . Top: retrieved transmission spectra for the atmosphere (Section 4.4.1) and starspot (Section 4.4.2) scenarios for JWST NIRSpec G395H Visit 1, but now including the LDSS-3C observations from Diamond-Lowe et al. (2018). Both retrievals include a free offset between the ground-based and JWST data. Overplotted are the ground-based observations of GJ 1132 b with no offset applied (grey data points), including the median retrieved offset from the starspot scenario (341 ppm; orange data points), and including the median retrieved offset from the atmosphere scenario (155 ppm; green data points). Bottom: posterior distributions including the ground-based data (green/orange for the atmosphere/starspot scenarios) compared to the results from Figure 4 with only JWST data (purple/blue for the atmosphere/starspot scenarios). The retrieval results minimally change when adding the ground-based observations.  \n<!-- image -->'}