The Featureless HST/WFC3 Transmission Spectrum of the Rocky Exoplanet GJ 1132b: No Evidence for a Cloud-free Primordial Atmosphere and Constraints on Starspot Contamination

Orbiting an M dwarf 12 pc away, the transiting exoplanet GJ 1132b is a prime target for transmission spectroscopy. With a mass of 1.7 M ⊕ and radius of 1.1 R ⊕, GJ 1132b’s bulk density indicates that this planet is rocky. Yet with an equilibrium temperature of 580 K, GJ 1132b may still retain some semblance of an atmosphere. Understanding whether this atmosphere exists and its composition will be vital for understanding how the atmospheres of terrestrial planets orbiting M dwarfs evolve. We observe five transits of GJ 1132b with the Wide Field Camera 3 (WFC3) on the Hubble Space Telescope (HST). We find a featureless transmission spectrum from 1.1 to 1.7 μm, ruling out cloud-free atmospheres with metallicities <300× solar with >4.8σ confidence. We combine our WFC3 results with transit depths from TESS and archival broadband and spectroscopic observations to find a featureless spectrum across 0.7 to 4.5 μm. GJ 1132b therefore has a high mean molecular weight atmosphere, possesses a high-altitude aerosol layer, or has effectively no atmosphere. Higher-precision observations are required in order to differentiate between these possibilities. We explore the impact of hot and cold starspots on the observed transmission spectrum GJ 1132b, quantifying the amplitude of spot-induced transit depth features. Using a simple Poisson model, we estimate spot temperature contrasts, spot covering fractions, and spot sizes for GJ 1132. These limits, as well as the modeling framework, may be useful for future observations of GJ 1132b or other planets transiting similarly inactive M dwarfs.


INTRODUCTION
If our Solar System terrestrial planets and moons are any indicator, rocky exoplanets likely possess a diverse population of atmospheres.From the thick CO 2 -dominated atmosphere of Venus, out to the tenuous N 2 atmosphere of Pluto, the rocky worlds of our Solar System have undergone significant atmospheric evolution (for a review see Encrenaz & Coustenis 2018).Many terrestrial exoplanets have likely experienced similar transformations leading to a variety of atmospheres (e.g., Segura et al. 2003;Grenfell et al. 2007;Hu et al. 2012;Forget & Leconte 2014;Luger & Barnes 2015;Schaefer et al. 2016;Kite & Barnett 2020;Grenfell et al. 2020).It is also possible that some rocky exoplanets never lost their primary atmospheres.Instead, they continue to maintain a slight H/He envelope consisting of < 0.1% of their masses initially accreted from the planetary nebula (Owen et al. 2020).Characterizing these terrestrial atmospheres will be crucial for understanding the formation, evolution, and potential habitability of these rocky worlds.
Of interest for atmospheric characterization are the rocky planets orbiting nearby bright M dwarfs.
The smaller stellar radii of these stars translate directly into larger transit depths for Earth-sized planets passing in front of them (Nutzman & Charbonneau 2008;Winn 2010).Moreover, these systems tend to be compact, with the habitable zone generally within 0.2 AU of the star (Kopparapu et al. 2013).As the probability of a transit is inversely proportional to the semi-major axis, there is a higher likelihood of a M dwarf planet transiting than one orbiting FGK dwarfs.
One notable planet for atmospheric study is the Earth-sized rocky planet, GJ 1132b.Discovered with the MEarth survey (Berta-Thompson et al. 2015), this 1.13 ± 0.056 R ⊕ and 1.66 ± 0.23 M ⊕ rocky planet orbits a nearby inactive M3.5 dwarf (Hawley et al. 1996;Dittmann et al. 2017;Bonfils et al. 2018).GJ 1132b receives 19× more bolometric insolation than Earth, implying a global equilibrium temperature of about 580 K (assuming uniform heat redistribution and bond albedo of 0), making it far too hot to be habitable.GJ 1132b receives more stellar flux than Mercury, but its higher mass and surface gravity lead to the questions: Does GJ 1132b possess an atmosphere?And if so, what is its atmospheric composition?GJ 1132b was targeted by ground-based observations seeking answers to these questions (Southworth et al. 2017;Diamond-Lowe et al. 2018).In broadband transit depths, Southworth et al. (2017) marginally detected hints of a low-metallicity atmosphere.This evidence was based primarily on a single anomalously deep z-band transit and required a 20% larger stellar radius than previously inferred.In contrast, Diamond-Lowe et al. (2018) observed a featureless spectrum between 700 and 1040 nm with an average precision of 100 ppm per wavelength bin (equivalent to 1.8 scale heights assuming a H/He mean molecular weight).Based on their higher precision results, they disfavor a cloud-free atmosphere of <10× Solar metallicity by volume.This suggests either a secondary atmosphere with a higher mean molecular weight, high-altitude aerosols, or no atmosphere around GJ 1132b.The broadband 4.5 µm Spitzer channel and the MEarth optical 0.7-1.0µm observations from Dittmann et al. (2017) do not by themselves strongly constrain possible atmospheres.
They do provide a light-curve-derived stellar density that confirms a stellar radius in agreement with Diamond-Lowe et al. (2018) and Berta-Thompson et al. (2015).Waalkes et al. (2019) searched for a deep Lyα ultraviolet transit of GJ 1132b using HST/STIS.
They detected no transits but place a 2σ upper limit of an exosphere at 7.3× the planetary radius.
From this, they argue that it is unlikely that GJ 1132b possesses an extended hydrogen envelope created from a leftover primary atmosphere or from the photodissociation of water.Based on the planet's total high-energy irradiation estimated from the stellar Lyα emission, Waalkes et al. (2019) calculate a mass-loss upper limit of neutral hydrogen to be 0.86 × 10 9 g/s, or one Earth ocean every 6 Myr.
Theoretical work shows that any atmosphere of a planet like GJ 1132b has been influenced by the intense X-ray and ultraviolet radiation and activity common to M dwarfs (France et al. 2016).
M dwarfs exhibit more XUV flux relative to bolometric than G dwarfs, making M dwarfs more efficient drivers of atmospheric escape (Lopez & Fortney 2013;Owen & Wu 2017).Analysis of UV observations of a similar planet host, LHS 3844, back these theoretical arguments while supporting the hypothesis that any atmosphere surrounding GJ 1132b is secondary in nature (Diamond-Lowe et al. 2021).Schaefer et al. (2016) modeled the evolution of GJ 1132b's atmosphere in order to predict potential atmospheric compositions, tracking both atmospheric escape and mantle outgassing.They find that GJ 1132b would have required at least 5% of its initial mass to be water in order to retain a waterdominated atmosphere today.As this would be unlikely for a planet forming interior to the snow line, Schaefer et al. (2016) argue that the most likely outcome is a tenuous atmosphere dominated by O 2 created by the photodissociation of H 2 O.However, most of the hydrogen and some of the oxygen atoms produced from this process are lost to space through hydrodynamic drag or absorption into a magma ocean.Schaefer et al. (2016) briefly comment on a CO 2 -dominated atmosphere.They note that adding significant CO 2 will prolong a surface magma ocean and enhance atmospheric loss.
However, future exploration and modeling is required to test this hypothesis.
Here we investigate the existence of an atmosphere around GJ 1132b by observing its transmission spectrum from the Hubble Space Telescope (HST) Wide-Field Camera 3 (WFC3) between 1.1 and 1.7 µm.These observations enable a sensitive search for water-vapor absorption, provide observational constraints on possible atmospheric compositions, and test whether or not GJ 1132b was able to retain a primordial hydrogen-rich atmosphere.
However, stellar activity, especially unocculted stellar spots and faculae, have the potential to either mimic or mask planetary water features in the WFC3 bandpass (Rackham et al. 2018;Zhang et al. 2018).GJ 1132, with a rotation period of 130 days, is not young (Newton et al. 2018) but does show variability due to rotating starspots, and Dittmann et al. (2017) attributed an offset between the MEarth optical and Spitzer infrared transit depths to unocculted stellar spots.We therefore also use available data to quantify the impact of GJ 1132's stellar activity on the interpretation of planetary transmission spectra observed by both HST /WFC3 and future instruments, such as the James Webb Space Telescope (JWST) or ground-based giant segmented-mirror telescopes (GSMTs).
We present this work in the following sections.We outline the observations in Section 2 and detail the analysis of the broadband and spectroscopic light curves in Section 3. We present an analysis of the TESS light curve of GJ 1132 in this section as well, providing additional constraints on the system properties and ruling out the transit of GJ 1132c.Section 4 discusses the transmission spectrum and the possible atmospheric compositions for GJ 1132b.We quantify stellar contamination in the transmission spectrum in Section 5 before concluding in Section 6.

Data
We observed five transits of GJ 1132b with WFC3/IR grism spectroscopy (GO #14758, PI Berta-Thompson).We used the G141 grism, providing stellar spectra at a resolution of R = λ/∆λ = 130 covering approximately 1.1−1.7 µm .Each visit consisted of one orbit used to settle WFC3's optical and detector systematics (which was discarded from analysis), one orbit before transit, one orbit in transit, and one orbit after transit.Phase constraints were set to maximize time in transit between Earth occultations, and the five scheduled visits together provide complete phase coverage of the transit, including both ingress and egress.
We gathered data with the 256x256 subarray, using the SAMP-SEQ=SPARS10 and NSAMP=15 settings for a total photon-counting exposure time of 88s per exposure.To minimize the overheads associated with detector readout and/or memory buffer dumps, we used a round-trip spatial scan in the cross-dispersion direction.A scan rate of 0.2./second kept the maximum fluence recorded on the detector a safe level of about 24,000 e-/pixel.For each exposure, 72s were lost to readout and scan resets, resulting in an overall photon-counting duty cycle of 55% for the 22 exposures gathered each orbit between Earth occultations.GJ 1132 is near the Galactic plane (b = 8 • ) and therefore in a crowded field.This presents a challenge for slitless spectroscopy due to the potential for multiple stars' spectra to overlap.In planning the observations, we picked a combination of ORIENT constraints and scan rate to prevent overlap of any significant background stars with GJ 1132's first-order spectrum during a single nondestructive read.

Reduction
We reduced each of the five WFC3 visits using the publicly available iraclis pipeline1 .We include a summary of the pipeline's steps here; Tsiaras et al. (2016a,b) provide a full description of the pipeline's methodology.Starting with the raw uncalibrated spatially scanned spectroscopic images, iraclis performs the following: a zero-read subtraction, non-linearity correction, dark correction, gain conversion, flat fielding, sky background subtraction, and bad-pixel and cosmic-ray interpolation.
iraclis calculates potential horizontal and vertical shifts in the scanned spectrum over time and determines the scan length for each image by approximating a Gaussian along the sum of the rows.
Based off the direct image of GJ 1132, wavelengths are assigned and the wavelength-dependent flat field is applied.iraclis reports times in units of HJD which we converted to BJD TDB using the time utilities code of Eastman et al. (2010) 2 .
We extracted the 1-D spectrum by applying a 166-pixel-wide aperture along the dispersion axis and summing along the cross-dispersion axis for each wavelength.As the average scan length for each image is about 170 pixels, we chose a slightly smaller aperture to minimize potential edge effects of the scan.However, slightly smaller and larger apertures (162 pixels and 172 pixels) had no notable improvements on the results.
Across the five visits, there is less than a 0.2 pixel shift in both the x− and y−directions of the scan.
A comparison of the residuals, created by subtracting off the best-fit models, to their corresponding x− and y− shifts demonstrated no correlation.We therefore did not apply any correction to these small shifts.We also found no correlation with either the scan length or background values to the residuals.Due to the crowded stellar field, we checked the direct and dispersed images of each visit for other nearby bright stars in the field of view.A comparison of the average spectrum between visits demonstrated less than a 1% flux variation across all wavelengths (Figure 1).We identified a weak 0.5% background source contaminating the GJ 1132 spectrum for each read, similarly noted in Mugnai et al. (2021).This background source contamination imparts a 12 ppm systematic suppression of the overall transit depth, smaller than the 34 ppm wavelength-binned spectroscopic transit depth uncertainties.
A comparison of each visit's absolute flux to the global median-combined total flux spectrum shows <0.2% deviation for visits 2−5 across 1.15 to 1.63 µm (Figure 2).Visit 1 demonstrates a <0.8% deviation with larger divergence at shorter wavelengths.This difference between the visits is likely a stellar effect.Visits 1 and 2 were observed 158 days apart (1.2× the stellar rotation rate of 130 days) while visits 2−5 were all observed within 64 days of one another (0.5× the stellar rotation).
To mitigate the slight flux variation across all visits, we opted to analyze each visit independently discussed in further detail below.

ANALYSIS
We analyzed the WFC3 data both as wavelength-integrated broadband light curves and as binned spectroscopic light curves.We also analyzed archival TESS data of GJ 1132, enabling both a new optical transit depth for GJ 1132b and a search for transits of the radial-velocity detected planet GJ 1132c (Bonfils et al. 2018).

HST Broadband Light Curve Analysis
From the iraclis-reduced spectra, we created broadband light curves for each visit by summing the flux from 1.15 to 1.63 µm.Uncertainties are a combination of the photon shot noise from the star, sky background, dark current, and read noises added in quadrature.As observed in multiple WFC3 light curves, each orbit demonstrated charge-trapping ramps, with the first orbit demonstrating significantly larger effects (e.g.Berta et al. 2012;Kreidberg et al. 2014;Wakeford et al. 2016;Zhou et al. 2017).We multiplied the physical charge-trapping instrumental model RECTE (Zhou et al. 2017) with the BATMAN transit model (Kreidberg 2015) and fitted each visit independently across all parameters (including both transit and instrumental systematic parameters).
For the RECTE portion of the model, we fitted for the starting number of slow and fast traps (E s,f ) as well as the number of traps between the orbits for slow and fast traps (∆E s,f ) recommended by Zhou et al. (2017).We also discovered that fitting for the total number of slow and fast traps (E s,tot and E f,tot ) significantly improved the fits for each visit.This is in contrast to Zhou et al. (2017), who recommended holding the total number of traps, trapping times, and trapping efficiencies constant.
We did find that varying trapping times and efficiencies did not improve our overall fits and held these constant to the values listed in Zhou et al. (2017).We included visit-long slopes and flux offsets to the charge-trapping model with these parameters dependent on the scanning direction.While Zhou et al.
(2017) note that one benefit of RECTE is the inclusion of the first orbit, we discovered that removing this orbit improved the overall fits for the three remaining orbits.Guo et al. (2020) also found this to be true in their WFC3 analysis of HD 97658b.RECTE is designed to approximate the charge trapping occurring per pixel.We therefore assumed an average per-pixel flux value by dividing the total flux in an image by the number of pixels contained in the aperture before fitting the systematics with

RECTE.
For the BATMAN transit portion of the model, we performed two independent versions of the transit fit: (1) varying a/R s and inclination and (2) holding these two parameters constant to the values determined by Dittmann et al. (2017) from high-cadence Spitzer data.For both versions, we fitted for the mid-transit time and R p /R s .Each visit only sampled the ingress or egress of the transit.We therefore held the quadratic limb-darkening coefficients constant, which were determined by LDTK (Parviainen & Aigrain 2015) in the WFC3 bandpass using GJ 1132's stellar parameters: log g = 5.049 (Dittmann et al. 2017), Z = −0.12(Berta-Thompson et al. 2015).We used the Gaia-determined stellar effective temperature of 3630 K (Gaia Collaboration et al. 2018) in calculating the limbdarkening coefficients.This T eff differs from the more accurate 3270 K Berta-Thompson et al. (2015) we adopt throughout the rest of this work; we re-ran one transit with the coefficients for this cooler temperature and found the parameter results for broadband and spectroscopic fits to be negligible.
The period remained fixed to that determined by Dittmann et al. (2017).We held the eccentricity constant at 0, which is supported by previous RV measurements from Berta-Thompson et al. (2015) and Bonfils et al. (2018).We also integrated the transit model over the WFC3 total exposure time of 103 seconds; accounting for this finite integration time is particularly important because it is comparable to the ingress/egress time (see Kipping 2010).
Our final model, combining both RECTE and BATMAN, for the flux F for,rev (t) as a function of both the forward and reverse scan direction and time t is therefore: (1) F for,rev and m for,rev represent the initial offset and visit long slope for the forward and reverse scanning directions respectively.Parameters denoted with {s, f } in the RECTE model refer to slow and fast trap parameters respectively, and parameters with 'ver' depend on the transit version in question.These parameters, combined with an error scaling parameter, lead to 15 and 13 fitted parameters for the first and second versions respectively.
We used the MCMC fitter emcee (Foreman-Mackey et al. 2013) to fit the above model to each of the five visits separately.An error scaling term (f scale ) was also included, which was multiplied to the flux uncertainties in order to allow the MCMC to inflate the uncertainties to achieve χ 2 r = 1 and propagate appropriately into the parameter distributions.Uniform priors were assumed for all parameters.We assigned 100 walkers with 15,000 steps.After an initial run, we check for any outliers that fell >5σ away from the median absolute deviation (MAD).These points were flagged and removed, and the fit was performed again.This accounted for two outlier points in visit 2 and one outlier point in visit 3. We binned the residuals from this fit and found a 1/ √ N Gaussian slope demonstrating that we properly isolated and removed systematic noise.
We decided to perform the two versions noted above as the transit duration of each visit is largely unconstrained.By varying a/R s and inclination for the first fit, we determined that all visits have less than a 2σ deviation from those values quoted in Dittmann et al. (2017).Therefore, we adopted Dittmann et al. (2017) values for these parameters (a/R s : 16.54 i: 88.68) for our second fit.By holding a/R s and inclination constant, we determined a more precise R p /R s and mid-transit time, though there is little deviation in the overall value between the two versions for both parameters.
Table 1 lists our best fit planetary parameters from both versions for the five HST visits.The second version parameters are used for all further analyses including the transmission spectrum analysis in Section 3.2.We included the parameters determined from the TESS light curve (Section 3.3) to Table 1 for comparison.The error scaling parameter (f scale ) fitted for during the MCMC routine and χ 2 r are listed as goodness-of-fit measurements.
We calculated the stellar density using the a/R s from each broadband transit and found that those of all but visit 2 were within 1σ of the 29.6±6.0g/cm 3 reported in Berta-Thompson et al. (2015).
Visit 2 has the least number of points defining the ingress or egress, making it a challenge to constrain the transit duration (and therefore stellar density).Southworth et al. (2017) measured a lower stellar density of 15.4 +4.8 −3.4 g/cm 3 , which gave them larger stellar and planetary radii; we find that all of our visits except for 1 and 2 have a >2σ deviation from this value.We therefore opt to use the stellar and planetary radii derived in Berta-Thompson et al. (2015) and Bonfils et al. (2018).
We plot the five broadband light curves with the best-fit models from the second version in Figure 3.
From top-down each panel shows: the normalized flux, the transit with the RECTE model removed, the systematics with the BATMAN transit model removed, and the residuals.Visit 1 has the largest deviation in flux between the forward and reverse scans, giving the ramp an up-and-down motion that is not as apparent in the other four visits.The origin of this difference is unclear, though the systematics model captures this effect.We do not plot the best fit from the first version, as the difference between the fits is slight and not apparent when plotted against each other.

HST Spectroscopic Light Curve Analysis
We created 22 spectroscopic light curves by summing the flux between 1.15 to 1.63 µm into 5-pixel bins (approximately 20 nm each).We analyzed these light curves using two different methods: fitting each light curve using RECTE and using the broadband divide-white method described in Kreidberg et al. (2014).held constant to that obtained from the second-version WFC3 broadband fit.This left R p /R s as the only free transit model parameter, along with the 10 systematic parameters from the RECTE model and an error scaling term.
Each spectroscopic light curve for each visit was modeled individually with the same emcee MCMC process as the broadband curves.The R p /R s for each visit and the corresponding wavelength bins are noted in Table 2. Inverse-variance weighted averages calculated from the five visits are listed in the last column of the table.Each spectroscopic light curve for all five visits is plotted in Figure 4 with the systematics-removed transit curve plotted on the left and the residuals for each bin to the right.We find little deviation in the RECTE model parameters from bin to bin, highlighting their mostly wavelength-independent nature.
Divide-White Method -We checked our RECTE analysis by applying the divide-white method to the spectroscopic light curves (Kreidberg et al. 2014).We created systematic residuals (third row in Figure 3) by dividing the broadband transit model from the broadband data.By assuming these systematic residuals are not wavelength-dependent, we divide these residuals from each spectroscopic light curve.We find that while the residuals captured the ramp-like effect of the data, they do not necessarily remove the difference in the forward and backward scanning offset.We therefore fit for two separate offsets, dependent on the scanning direction, as well as the transit depth and an errorscaling term.Again, we performed the same emcee MCMC routine as the broadband and RECTE analyses discussed above.
A comparison of the transit depths between the RECTE and divide-white methods demonstrated discrepancies that were <2σ for all but the second wavelength bins (Figure 5).The RECTE method always produced slightly larger error bars on the transit depths due to the larger number of free parameters.We opted to use the RECTE transit depths (R p /R s ) 2 for the rest of the analysis and discussion in this paper as this model is based on physical properties of the detector.

TESS Light Curve Analysis of GJ 1132b
GJ 1132 was observed by TESS during Sectors 9 and 10.In Sector 10, the star appeared on a noisier section of the CCD, demonstrating larger instrumental variation.For this analysis, we therefore opted to use only the Sector 9 data, including a total of 14 individual transits.We used the python package lightkurve (Lightkurve Collaboration et al. 2018) to download and analyze the two-minute-cadence observations.Starting from the PDC light curve, we detrended the flux by masking out all GJ 1132b transits and applying a median moving boxcar filter with a baseline of one day (32× the transit duration).As with the broadband WFC3 transits, we performed the same two versions of fits as above (varying a/R s and inclination and holding them constant) to all transits simultaneously.We include a baseline of 135 minutes (3× the transit duration) on either side of each transit.For both versions, we fitted for R p /R s , a transit ephemeris, quadratic limb-darkening coefficients, and a general offset.We assumed uniform priors on all parameters except for the limb darkening coefficients which were assigned Gaussian priors based on values calculated by LDTK in the TESS bandpass (Parviainen & Aigrain 2015).We include the error scaling parameter.Table 1 lists the best-fit parameters for both versions, and Figure 6 plots the best-fit model atop the folded  data.We found no significant difference between the two versions; the fitted values of a/R s and inclination were consistent with those from Spitzer (Dittmann et al. 2017).
We also analyzed individual TESS transits, again applying the two versions.Using the a/R s and inclination values for each transit, we calculated the transit duration for each epoch.We see no evidence for transit duration variations or transit depth variations in the TESS Sector 9 data.A brief inspection of Sector 10 showed large variations in each due to significant instrumental noise in the data even after the median detrending, supporting our decision to remove this sector from our analysis.
We measured the Sector 9 mid-transit times for each transit, holding a/R s and inclination constant.
We combined these times to the broadband mid-transit times for each HST visit as well as the published times from MEarth (Berta-Thompson et al. 2015;Dittmann et al. 2017), Spitzer (Dittmann et al. 2017), and MPG telescope (Southworth et al. 2017).Diamond-Lowe et al. ( 2018) noted that their mid-transit times for GJ 1132b were consistently 2 minutes early for all four of their observed  mid-transit times for each epoch, we fit a straight line through these points providing an updated period of 1.6289299 ± 5x10 −7 days and a transit ephemeris of 2457184.55747± 0.00012 BJD TDB .
Table 3 reports the mid-transit times for HST broadband curves with the inflated uncertainties and times from Sector 9 of TESS.A search for transit timing variations (TTVs) yielded a non-detection with the residuals matching a flat line with χ 2 r = 1.17 for 53 degrees of freedom (Figure 7).

TESS Light Curve Analysis of GJ 1132c
Discovered with radial velocities by Bonfils et al. (2018), GJ 1132c has not yet been completely vetted for the possibility of transits.Bonfils et al. (2018) estimate a 1% transit probability.No obvious transits were seen in MEarth or Spitzer data, but they lacked complete and continuous coverage of the possible times of transit.We use TESS's continual coverage to search for the possibility of a GJ 1132c transit.We predict two possible transit events of GJ 1132c to occur in each Sector.As with  our GJ 1132b TESS analysis, we limited our transit search to Sector 9. We masked out all GJ 1132b transits in the median boxcar detrended data.We assumed Gaussian priors on the a/R  time and period parameters using values determined by Bonfils et al. (2018) and set uniform priors on R p /R s and the inclination.We allow for grazing transits by letting the impact parameter to vary from 0 to 1+R p /R s .We also allow for the possibility of negative transit depths.The quadratic limb darkening coefficients were held constant to those values determined from the GJ 1132b TESS fits.
We performed an MCMC fit using 100 walkers and 15,000 steps with 5,000 of those removed for burn-in.We determined an R p /R s of 0.024 +0.011 −0.046 , a result consistent with zero.A visual inspection of the data shows no indication of a obvious transit (Figure 8).From the posterior, we place a 3σ upper limit on R p /R s at 0.081, marginalized over all other parameters.This corresponds to a planet radius of 1.84 R ⊕ .With a minimum mass of 2.64 M ⊕ (Bonfils et al. 2018), it is possible that this planet possesses a radius smaller than what we can detect with two TESS transits.More data from the TESS extended mission will be required to completely rule out a GJ 1132c transit.We also

O -C
Figure 8.The TESS Sector 9 light curve of GJ 1132, folded assuming the best-fit period (8.922 ± 0.009 days) from our TESS analysis of GJ 1132c, with the best-fit maximum likelihood model plotted in light blue.We find no significant evidence that GJ 1132c transits.
add an additional caveat that a box least squares (BLS) algorithm would provide a more rigorous search of the data, especially as the RV parameters on the period and expected mid-transit time are not well-constrained.However, we calculate that from GJ 1132c's a/R s (Bonfils et al. 2018), its orbit would need to be inclined to at least 89.88 degrees (1.2 degrees from GJ 1132b's orbit) in order to transit.While this is close to co-planar with its transiting neighbor, the tight constraint on GJ 1132c's inclination for transits to occur makes it more likely that it is inclined such that it does not transit.spread, except for wavelengths 1.25 and 1.40 µm.Removing the outlying visits in these bins had no appreciable effect on the final spectrum.
Before combining the five visits into a weighted-average spectrum, we checked that the fitted transit depths in each bin are drawn from the same distribution and there are no significant outliers.We determined that every individual transit depth fell within 2σ of the weighted-average depth in their respective bin, except for the 1.25 and 1.40 µm bins (Figure 9).Both bins showed a 2.2σ spread with visits 3 and 4 yielding the outlying points respectively.A comparison between the average transit depths with and without these visits yielded no difference on our final results.We therefore keep all five visits when calculating the average transit depth for each wavelength bin.and log(g): 5.0 and translated to the MEarth bandpass using the same method described in Irwin et al. (2018).Parameter uncertainties were scaled to account for overdispersion globally rather than individually per night.This re-analysis yielded a deeper transit depth of 0.00246 ± 0.00011, a 3σ difference than the shallower transit found in Dittmann et al. (2017).It is also statistically the same (<1σ) as both the Spitzer 4.5 µm and the WFC3 depths.As such, this updated depth indicates no significant stellar contamination between the shorter optical and longer infrared wavelengths.
•  2018) is deeper than any of their spectroscopic depths while also being statistically similar to the depths derived from the WFC3 observations (Figure 11).We hypothesize that this difference is due to instrumental and systematic effects and not astrophysical in nature.We add a fitted offset of 0.000205 to the spectroscopic transit depths such that they best fit the various metallicity models which we discuss below.
• Swain et al. ( 2021) used these HST/WFC3 data to infer a Rayleigh scattering slope and a possible HCN and CH 4 absorption feature.We do not observe this in our analysis, and we use our featureless WFC3 spectrum for the comparison with the other data sets.
• Mugnai et al. (2021) also used these HST/WFC3 data but determined a featureless transmission spectrum in agreement with our transit depths.Again, we opt to use our WFC3 depths for the rest of the discussion.

Cloud-Free Atmospheric Models
We present the WFC3 weighted-average transit depths for GJ 1132b in Figure 12.The planet's transmission spectrum is best explained by a flat or featureless spectrum with a χ 2 r of 1.02 with 21 degrees of freedom.The average transit depth uncertainty is 34 ppm per 20 nm wavelength bin.On GJ 1132b, an H/He-dominated atmosphere with a mean molecular weight of µ = 2.2 would have a scale height of 172 km; the WFC3 transit depth uncertainties correspond to an effective resolution of 0.29 of those scale heights.In turn, this precision corresponds to 1 scale height for an atmosphere with µ = 7.6.We expect absorption features in the atmosphere to possess sizes of several scale heights.From this precision and the featureless spectrum, we can confidently rule out atmospheres with mean molecular weights less than µ = 7.6.
For interpreting these transit depths, we generated model transmission spectra using Exo-Transmit  Figure 12.The WFC3 transmission spectrum of GJ 1132b, averaged over the five visits of WFC3, plotted against cloud-free atmospheric models of varying metallicities for comparison.For each model, we provided the χ 2 r and the p-value evidence against the model.We translated the p-value to the number of σ (nσ) confidence that we can reject this model.The spectrum is most consistent with a flat line indicating a high mean molecular weight atmosphere or no atmosphere.We could easily detect an atmosphere with a metallicity less than 300× Solar by volume (or mean molecular weights less than 8.9 amu).
As the 1 bar pressure radius is unknown, we use GJ 1132b's planet radius to generate Exo-Transmit models for each metallicity before re-scaling them by a multiplicative factor less than 1 in order to best fit the HST transmission spectrum (Kempton et al. 2017;Diamond-Lowe et al. 2018).
From the WFC3 observations, we confidently rule out atmospheres with metallicities <300× Solar by volume at >4.8σ significance.This corresponds an atmospheric mean molecular weight >8.9 amu (if an atmospheric exists).For comparison, a 1× Solar atmosphere has a mean molecular weight of  12 except with 100% molecular composition atmospheric models for comparison.
We are unable to rule out any of these atmospheric models as their higher mean molecular weights produce smaller transit depth variations than the precision of the data points.
atmospheres of 600× Solar and 1000× Solar metallicities fit the WFC3 spectrum with χ 2 r of 2.6 (p-value of 0.0025) and 1.3 (p-value of 0.14) respectively.
However, given GJ 1132b's equilibrium temperature of 580 K (Bonfils et al. 2018), it is more likely that this planet possesses either a CO 2 -dominated atmosphere similar to that of Venus, a tenuous O 2 atmosphere as discussed in Schaefer et al. (2016), or no atmosphere at all.When comparing atmospheric models of 100% H 2 O, O 2 , and CO 2 with the WFC3 data, we calculate χ 2 r of 1.8 (p-value of 0.01), 1.1 (p-value of 0.34), and 1.2 (p-value of 0.26) respectively.From these statistics and as illustrated in Figure 13, we cannot rule out any of these scenarios with the current WFC3 precision.
We note, however, that GJ 1132b's equilibrium temperature makes it very unlikely that this planet possesses a H 2 O-rich atmosphere.

Models With Clouds and Hazes
Featureless exoplanet transmission spectra can be attributed to two different causes: a higher mean molecular weight atmosphere (and thus smaller scale height) or a high-altitude aerosol layer muting or blocking absorption features.While GJ 1132b's temperature and terrestrial nature suggest a higher mean molecular weight atmosphere, the WFC3 observations alone cannot rule out the possibility of a H/He dominated atmosphere combined with condensate clouds such as ZnS or KCl (Morley et al. 2013) or organic photochemical hazes.We therefore determine the altitude required for an aerosol layer to mute all features in the WFC3 transmission spectrum.Assuming a 1× Solar composition, we place aerosol layers at increasing altitudes until we obtain a spectrum that agrees with the data (p-value of 0.001; process described in detail in Libby-Roberts et al. 2020).We find that this layer must lie at pressure levels <0.4 mbars in order to mute the features of a H/He dominated atmosphere.
For comparison, Kreidberg et al. (2014) determine the cloud layer on the sub-Neptune GJ 1214b lies at pressure levels <0.01 mbars assuming a low mean molecular weight atmosphere.
We test the possible existence of an aerosol layer above a 1× Solar atmosphere by calculating GJ 1132b's mass loss rate for a H/He primordial envelope.Owen et al. (2020) show that low-mass planets similar to GJ 1132b can accrete around 1% of their masses in H/He gas during formation.With a current mass of 1.66 M ⊕ , this would translate to a starting H/He envelope mass of 0.0166 M ⊕ .
Waalkes et al. ( 2019) approximated a neutral hydrogen mass-loss rate of 0.86 ×10 9 g/s.Assuming no secondary sources of hydrogen and a constant mass-loss rate, GJ 1132b would lose such an atmosphere in 3.7 Gyr.Berta-Thompson et al. (2015) approximate this system to be >5 Gyr.Moreover, the stellar flux from GJ 1132 should have been higher in the past which would have increased this rate.
Therefore, it is very unlikely that GJ 1132b could possess a present-day cloudy primordial H/He-rich atmosphere.

GJ 1132b's Atmospheric Composition Determined From Archival Data
A direct comparison between all archival data to the various models shows that every atmospheric model (featureless or otherwise) are ruled out with > 5σ confidence (Figure 14).This is largely driven by the offset of Diamond-Lowe et al. ( 2018) spectroscopic points from the other archival transit depths.Given the discrepancy between the spectroscopic and broadband LDSS3C transit depths, we treat the relative changes in the transit depths within the spectrum as accurate, but not as the absolute depth.We vertically shift these points such that they best fit the different atmospheric models along with the TESS, MEarth, and Spitzer broadband points and the WFC3 spectrum (Figure 15).We recalculate the goodness-of-fits for all atmospheric models after fitting for this offset.The addition of literature data across 0.7 − 4.5 µm enables us to rule out 1x, 100x, and 300x Solar metallicities by volume with confidences of 17.6σ, 10.3σ, and 5.6σ respectively.A flat line remains the best-fit model with a χ 2 r of 1.01 (p-value: 0.45), while 600x and 1000x Solar metallicity models have χ 2 r of 2.0 (p-value: 0.00035) and 1.5 (p-value: 0.033) respectively.Notably, the Spitzer 4.5 µm point disagrees with the 600x Solar metallicity model with >2σ.However, it remains the only observation at wavelengths longer than the 1.63 µm cutoff from the WFC3 spectrum.Future JWST observations will be able to fill in this wavelength gap with a precise spectrum.

STARSPOT CONTAMINATION
Spatial inhomogeneities on the surface of a star can influence planetary transit depths in a wavelength-dependent fashion, potentially mimicking the signal of absorption through a planetary atmosphere.A planet transiting across a dark spot or a bright facula can show a wavelength-dependent bump (Pont et al. 2008;Espinoza et al. 2019) that could modify transit depths if not identified and corrected.More perniciously, even a planet transiting an unspotted patch of photosphere will be blocking starlight that is unrepresentative of the average for the star as a whole.This introduces features in observed transit depths (and transmission spectrum) that express the difference between the spotted and unspotted surface (Czesla et al. 2009;Sing et al. 2011).Termed the "transit light source effect" (TLSE) by Rackham et al. (2018), this unocculted starspot phenomenon poses a particular challenge for infrared observations of planets transiting cool stellar hosts.M dwarfs have strong molecular features in their photospheres that vary in-and-out of spots, meaning they can introduce spurious transit depth features from molecules like water that might also be expected in a planet's transmission spectrum (see also Zhang et al. 2018;Wakeford et al. 2019).Here we apply a simple model of the TLSE to the flat WFC3 transmission spectrum in order to infer the starspot properties of the host star GJ 1132 and quantify the impact of the TLSE on current and future observations of 1X Solar ( 2 : 11.1; pval: 1.7e-69; n : 17.6) 100X Solar ( 2 : 5.6; pval: 1.3e-27; n : 10.8) 300X Solar ( 2 : 3.5; pval: 5.5e-13; n : 7.1) 600X Solar ( 2 : 2.9; pval: 1.3e-9; n : 5.9) 1000X Solar ( 2 : 2.7; pval: 4.7e-8; n : 5.3) Flat Line ( 2 : 2.9; pval: 2.5e-4; n : 5.9) Goodness-of-fits for each model were calculated using the LDSS3C spectrum, TESS, MEarth, and Spitzer photometry, and the WFC3 depths from this work.The noisy z−band flux point (Southworth et al. 2017) was cropped out to ease model visualization on the y−axis.From the χ 2 r calculated with as-published depths, it appears that no planetary atmosphere model explains all transit depths; however, this is largely driven by the disagreement in the optical LDSS3C spectroscopic depths.

Libby
GJ 1132b or other similar planets.We outline the modeling framework, enumerate the inputs and assumptions we adopt, and present TLSE constraints.

Starspot Model
The wavelength-dependent transit depth of a planet transiting a spotless star can be written as , where the first term represents the transit depth of the planet at a reference radius (such as a rocky surface or a 1-bar level in its atmosphere) and ∆D(λ  larger by n(λ) scale heights at different wavelengths.In this section, we instead aim to model D(λ) as an atmosphere-less planet transiting a spotted star.We rewrite the wavelength-dependent depth as D(λ) = (R p /R s ) 2 + ∆D(λ) spot , where the spot-induced depth variations (∆D(λ) spot ) can be either positive or negative.Approximately following the derivations outlined in Rackham et al. (2018) and Zhang et al. (2018), the spot-induced depth variation can be estimated as with effectively three unknowns: the flux ratio S(λ) spot /S(λ) unspot , the global spot covering fraction f , and the transit-chord spot covering fraction f tra , which are summarized in Figure 16 and below.
• The spectral fluxes S(λ) spot and S(λ) unspot of the spotted and unspotted photospheres, respectively, represent the flux (in W/m 2 /nm or related units) emitted by two different types of stellar photosphere.For this work, we use solar-metallicity model stellar spectra from Husser et al. ( 2013)3 for a surface gravity of log g = 5.0, which closely matches that of GJ 1132, so the temperature (T spot or T unspot ) uniquely determines the spectrum.To maintain generality for modeling spots or faculae, T spot is allowed to be either cooler or warmer than T unspot .By using fluxes instead of angle-dependent intensities, we are ignoring the effects of limb-darkening but still sensitive to the effect of starspots averaged over the entire stellar disk.
• The average spot covering fraction f represents the long-term average of the fraction of the Earth-facing stellar disk covered with spots.The instantaneous spot covering fraction facing Earth at any particular time t can be written as f ±∆f (t), where ∆f (t) represents the variability due to rotation (or, potentially, spot evolution).For example, a star with an average spot covering fraction of f = 20% might exhibit an Earth-facing spot covering fraction ranging from 19% to 21% as it rotates through one or more periods.∆f (t) therefore has a semi-amplitude of 1%.It is important to separate the time-variable term from the average spot covering fraction, as the TLSE depends on the total spot covering fraction in and out of the transit chord, not just the (always smaller) fraction producing out-of-transit modulations.The model for ∆D(λ) spot described here represents the time-averaged impact of starspots on the transit depth rather than the variability from epoch to epoch.This effectively assumes that transits will be observed over multiple epochs and/or that the static signal will dominate over the time-variable contribution to the transit depth.
• The transit chord spot covering fraction f tra represents a similar time average along the narrow slice of the stellar disk transited by the planet; the amplitude of ∆D(λ) spot depends on the difference in the spot covering fraction between the disk and the surface actually blocked by the planet.For this work, we assume that f tra = 0. Motivated by the lack of obvious spot-crossing = the spectrum of the spotted photosphere = the spectrum of the "unspotted" photosphere The average spot-covering fraction = = the spot-covering fraction facing toward Earth, which varies as the star rotates The average spot-covering fraction along the transit chord = Figure 16.Definition of the variables used for modeling the impact of stellar surface inhomogeneities on the inferred transmission spectrum.In this framework, 'spots' can represent patches that are either cooler or hotter than the surrounding stellar photosphere.
features observed in transits, this assumption considerably simplifies Equation 2. However, the f we infer through modeling will more accurately represent how spotted the average disk surface is relative to the transit chord.It is effectively assuming that the stellar flux along the transit chord can be well represented by a single-temperature PHOENIX atmosphere; the rest of the stellar disk comprises a f -weighted average of spectra at two different temperatures.

Starspot Data and Priors
To connect this conceptual framework directly to the specific case of GJ 1132b transiting its midto-late M dwarf, we apply three data constraints and three priors.
Data Constraint (1) = The semi-amplitude of out-of-transit modulations -In this model, as a star rotates, the time-variable integrated stellar flux S(λ) + ∆S(λ, t) should change with time according to which when integrated over any particular bandpass can be compared directly to the semi-amplitude of observed photometric modulations.The star GJ 1132 has been monitored nearly continuously by the MEarth Observatory (Newton et al. 2018) since before the discovery of GJ 1132b.These data trace the photometric variability ∆S(λ, t)/S(λ) of the star due to rotation or starspot evolution, integrated over the MEarth RG715 filter bandpass (roughly 700−1000 nm).MEarth out-of-transit monitoring data for GJ 1132, spanning February 2014 through August 2019, are shown in Figure 17 4 .For the data shown here, all high-cadence transit follow-up observations were excluded.A correction for time-variable precipitable water absorption in Earth's atmosphere was applied using a linear decorrelation against the MEarth common mode (see Berta et al. 2012;Newton et al. 2016).
Some residual trends are not perfectly corrected through this process, as seen by the disagreements where data from two MEarth telescopes overlap.From the subset of these data that were available at the time, Berta-Thompson et al. (2015)  The more recent data (Figure 17) continue to support these inferences of a stellar rotation period of roughly 120-130 days.The variation in the light curve's appearance from season to season over 5 years of observations clearly indicates significant starspot evolution and/or differential rotation.We estimate by eye the semi-amplitude of ∆S(λ, t)/S(λ) to be about 1.0±0.2% in the MEarth bandpass.
Though TESS and Spitzer also provide precise out-of-transit monitoring, the datasets span durations much shorter than the star's rotation period and therefore cannot be used to estimate the full range of rotational variability of the star.The out-of-transit spectra from this WFC3 program could be used to spectroscopically constrain ∆S(λ, t)/S(λ) at 1.1-1.7 µm, but Figure 17 shows the five WFC3 epochs do not sample the star at its full range of variability, so we exclude them from the analysis.with σ as the Stefan-Boltzmann constant.Without this, arbitrary combinations of spot covering fractions and temperatures could produce stars with extremely unrealistic integrated surface fluxes.
Data Constraint (3) = The transit depths -In this model, D(λ) can be calculated with Equation 2and compared directly to observed transit depth at any wavelengths.For this work, the main transit depths we include are the measured WFC3/G141 transit depths reported in Table 2.We also include the spectroscopic LDSS3C depth measurements from Diamond-Lowe et al. (2018) because they are the most constraining optical measurements available.As in Section 4.1.3,we treat the absolute depth of the LDSS3C as unknown.When comparing to starspot models, we include a multiplicative offset allowing these depths to shift up and down together as a free parameter.We marginalize over this offset in all our starspot inferences.These data are shown in Figure 18, with the maximum-likelihood offset applied to the LDSS3C depths.
Prior (1) = The connection between ∆f (t) and f -In this model, the total spot covering fraction f and the variability in the Earth-facing fraction due to rotation ∆f (t) are entirely decoupled.However, through a suite of geometric simulations, Rackham et al. (2018Rackham et al. ( , 2019) ) determined the amplitude of photometric rotational modulations generally scales as f 0.5 .Such a scaling emerges naturally by treating ∆f (t) as governed by a Poisson process set by the average number of spots covering the stellar disk.If we define f 1 as the average spot covering fraction of one individual spot, the number of spots over the entire stellar disk would be N = f /f 1 , and the expected random variation in the number of spots visible at any one time would be √ N .This implies the expectation value for the level of rotation variability ∆f (t) should be This broadly encapsulates the behavior that high-amplitude rotation variability can emerge from either larger spots (higher f 1 at fixed f ) or a more thoroughly spotted star (higher f at fixed f 1 ).To capture this relationship, we apply a prior that x = [f − ∆f (t)]/f 1 = N − ∆f (t)/f 1 should be drawn from a Poisson probability distribution with N as the expectation value, written as  2018) assume quasi-isotropic symmetries to estimate rotational modulations scale as f 0.5 .The typical spot size distribution cannot be determined robustly if spots instead congregate to fill active latitudes or polar caps (Guo et al. 2018).
Prior (2) = Uniform logarithmic priors on f , ∆f (t), f 1 -As the values of the spot covering fraction parameters could plausibly span orders of magnitudes, we wish for their priors to be uniform for the logarithm of the value.Therefore, we apply a logarithmic prior of the form P (x) ∝ 1/x on each of f , ∆f (t), and f 1 .With this prior, the probability will be uniform across logarithmic intervals, meaning the prior probability of a parameter falling in the range 10 −3 -10 −2 is the same as for the range 10 −2 -10 −1 .Furthermore, as the effect of covering f > 0.5 of the star with cool spots can be reproduced with f < 0.5 coverage by hot spots, we place a strict prior of f < 0.5 to avoid those (effectively duplicate) solutions.We also apply a prior requiring f 1 < ∆f (t) < f .
Prior 3 = The range of allowable spot temperatures -In this model, spot temperatures could be arbitrarily higher or lower than the unspotted surface.However, to avoid introducing an artificial bias toward hot spots due to the practical limitation that the PHOENIX spectra we use are not available below 2300 K, we impose a symmetrical prior that 0.75 < T spot /T eff < 1.25.This range does remove some otherwise viable models, but it still includes most of the spot temperature contrast ratios inferred for cool dwarfs in the Berdyugina ( 2005

Starspot Inferences
We sample from the posterior distribution for the parameters T unspot , T spot , f , ∆f (t), f 1 , and R p /R s (as well as a nuisance offset for the LDSS3C depths).We use the emcee (Foreman-Mackey et al.

2013
) ensemble sampler with 100 walkers.We sample for a total of 20,000 steps and exclude the first 10,000 steps as burn-in.Visual inspection showed that parameter distributions already appear converged by 1000 steps.
For a first experiment, we assume GJ 1132b's atmosphere has such a high mean molecular weight that its transmission spectrum can be entirely neglected.Thus the planet's radius is effectively constant across wavelength and all transit depth variations are due to the effect of unocculted starspots.
Figure 18 visualizes the resulting constraints on GJ 1132's starspot properties.We plot the combinations of fitted parameters that most clearly explain phenomena in the model, and color each sample by the value of its starspot temperature ratio T spot /T unspot so that hotter and colder spots can be traced throughout.We summarize some key conclusions from Figure 18 and this analysis as follows.
• The data for GJ 1132 can be closely matched by the starspot model outlined above.At the peak of the posterior distribution, the data yield χ 2 = 46.4 for 41 data points (1 MEarth semi-amplitude, 1 effective temperature, and 17 LDSS3C and 22 WFC3 depths) and 7 free parameters.
• The existing data are insufficient to distinguish whether GJ 1132's most important unocculted surface features (either spots of faculae/plage) are cooler or warmer than the surrounding photosphere.More precise optical data, more reliable knowledge of the absolute depths offset between LDSS3C and WFC3, and/or precise spectroscopic measurements of the rotational modulation semi-amplitude could distinguish between cold or hot spots.
• The distribution of ∆f (t) and T spot /T unspot is shaped by the need to match the semi-amplitude of the rotational modulations; MEarth's observed 1% variability can be matched either with stronger spot contrasts or more dramatic variations in the spot covering fraction.Only hot spots can extend to ∆f (t) < 1% and still match the 1% photometric variability because only they can produce spot flux contrasts > 100%, whereas the most that a completely dark cold spot could change a region's surface flux by would be 100%.
• The distribution of ∆f (t) and f is shaped additionally by the flat transit depths.f would otherwise uniformly fill the upper left triangle above f = ∆f (t), but the lack of observed spectral features in the transit depths rules out models with low variability ∆f (t) (which also have strong temperature contrast) and high average covering fraction f .
• The average number of visible spots N = f /f 1 is mostly < 10, as shown by dashed lines plotted on the distributions for f and N .Models with large numbers of very small spots would require too high of f for a given ∆f (t) and thus introduce strong spectral features into the transit depths.The values of f = 0.044 +0.094 −0.027 are generally 2 − 3× the values of ∆f (t) = 0.018 +0.029 −0.008 .
The typical single-spot covering fraction f 1 = 0.010 +0.016 −0.006 ; as a fraction of the 2π steradians of the Earth-facing side of the star, this would correspond to typical spot radii of 1.1 +1.3 −0.8 degrees.
These sizes are about 2× larger than the Sun's largest spots at solar maximum (Mandal et al. 2017) and similar to the largest active regions recorded on the Sun (Hoge 1947;Nicholson 1948).If the planet transited one of these spots at disk center, the crossing event would last about a minute and likely fall entirely within a single WFC3 exposure.We did observe a few single-point outliers in the broadband transit light curves (Figure 3), but as there were also outliers in the out-of-transit baseline we cannot uniquely attribute them to starspot crossing events.
• The true planet-to-star radius ratio R p /R s changes depending on whether the unocculted spots are cold or hot.If the spots are cold, the transit chord is brighter than the rest of the star, so the transit depth over-estimates R p /R s (and vice versa).Marginalizing over both cold and hot spots, we find the true R p /R s = 0.0490 +0.0005 −0.0003 .For GJ 1132b, the flat transit spectrum constrains the spot-induced uncertainties on R p /R s to be broadly similar to the statistical uncertainties on this quantity, and negligible compared to the uncertainties on the stellar radius.
The above conclusions connecting rotational variability to the average spot covering fraction rely on an assumption that the longitudinal distribution of spots is random and isotropic.This could be broken by, for example, a large polar spot or an extremely dense zonal band.Such extreme symmetries would exhibit little photometric variability but imbue conspicuous spectral features into the WFC3 depths.We see the opposite (strong photometric variability and flat WFC3 depths), suggesting the Poisson-based approach used here may be valid.
For an additional experiment, we ask whether an unfortunate arrangement of starspots could artificially erase the transmission spectrum of a H/He-rich atmosphere by producing a starspot signal that exactly cancels the extra transit depth from the planet's atmosphere.To answer this question, we produce a fake dataset into which we inject a flipped version of the model transit depths from the above 100× Solar metallicity atmosphere.If the starspot model could explain those flipped depths, the flat observed spectrum could just be a coincidence.We repeat the sampling procedure as above but for simplicity exclude the LDSS3C depths.In this fit, hot spots with strong temperature contrasts and large covering fractions can produce inverted water features that qualitatively mimic those needed to mask the planet's transmission spectrum (as also seen in TRAPPIST-1, Zhang et al. 2018).
However, at the peak of the posterior, χ 2 = 59.2 for 24 data points (1 MEarth semi-amplitude, 1 effective temperature, and 22 WFC3 depths) and 6 free parameters.This is a poor match to the data (p-value of 3 × 10 −6 ), mostly driven by the mismatch in shape between the cool water features in the transmission model and the hot water features in the stellar spectra.From this test, we disfavor a scenario where GJ 1132b has a clear H/He-dominated atmosphere that happens to be masked by starspot contamination.
Another related way of partially cancelling out strong transmission features is the spectral resolution-linked bias (Deming & Sheppard 2017).Transit depths observed at low spectral resolution will be biased toward wavelengths where the star is brightest within each bin.If strong planetary absorption features align in wavelength with deep stellar absorption features, as water lines would do for GJ 1132b and GJ 1132, transmission features will be suppressed.Deming & Sheppard (2017) show the feature suppression to be at the level 12% for TRAPPIST-1 (500 K cooler than GJ 1132) at WFC3 wavelengths; we do not account for this effect anywhere in our analyses.
Currently available data suggest that unocculted starspots affect GJ 1132b's transit depths at about the 100 ppm level in the optical and about the 10 ppm level for wavelengths longer than 1 µm.
Starspots do not pose a serious problem for the HST/WFC3 transmission spectrum of this target, but future analyses of transit observations with JWST or other large telescopes should be aware of the potential contaminating influence of the star's mottled surface.

CONCLUSION
We investigated the HST/WFC3 transmission spectrum of GJ 1132b, a rocky super-Earth orbiting a nearby bright M dwarf, over five separate visits.We determined a featureless spectrum to a precision of 34 ppm in 20 nm wavelength bins spanning 1.15-1.63µm.From this result, we rule out the presence of a cloud-free H/He-dominated atmosphere with a mean molecular weight less than 8.9 amu at 4.8σ confidence.High-altitude aerosol layers at pressures less than 0.4 mbars could potentially flatten the features of a solar composition atmosphere.Using the predicted mass-loss rate estimated by Waalkes et al. (2019), we find that a primordial H/He atmosphere composing 1% of GJ 1132b's mass could be removed in 3.7 Gyrs by the current high-energy irradiation from the star GJ 1132.
With an age >5 Gyrs and the much more intense radiation environment it experienced in the past, GJ 1132b is unlikely to possess any primordial H/He dominated atmosphere.We therefore conclude that the most likely scenario for GJ 1132b is that it has a high-metallicity secondary atmosphere, such as Venus' CO 2 -dominated atmosphere, the eroded O 2 -dominated atmospheres predicted by Schaefer et al. (2016), or no atmosphere at all.Current data cannot speak to the total mass in such a secondary atmosphere, which could range from a Venusian-like thick atmosphere to Mercurial-like tenuous exosphere.This conclusion is contrary to the recent work by Swain et al. (2021), who use a different analysis of the same data set to infer an H/He-rich atmosphere with a Rayleigh scattering slope and HCN and CH 4 absorption.Our results are consistent with the featureless spectrum found by Mugnai et al. (2021).
We analyzed the GJ 1132 TESS light curve and determined a broadband transit-depth and orbital parameters for GJ 1132b consistent with those from the WFC3 transits and previous Spitzer measurements in Dittmann et al. (2017).A search for GJ 1132c transits in the TESS data was unsuccessful, yielding a 3σ upper limit of 1.84 R ⊕ on the radius of the planet if it transits.Indeed, if GJ 1132c was completely co-planar with GJ 1132b, it would not transit.Instead GJ 1132c requires an orbital inclination of 89.88 degrees (1.2 degrees difference from GJ 1132b) or more for a transit to occur along our line-of-sight.Bonfils et al. (2018) note that there is a less than 1% chance this planet transits, supporting our non-detection in the TESS light curve.
Combining our GJ 1132b WFC3 spectrum with the TESS broadband depth, our updated MEarth depth, and other archival GJ 1132b transit depths (Dittmann et al. 2017;Diamond-Lowe et al. 2018) yielded a transmission spectrum covering 0.7-4.5 µm for this planet.We determined that the The best-fit model across all points remained a featureless flat line with a χ 2 r of 1.01 (p-value: 0.45).
Future JWST observations or ground-based observations from the upcoming ELTs will be helpful in discerning the existence of an atmosphere around GJ 1132b and, if one exists, determining its composition.
We explored the influence of unocculted spots on the measured transit depths to assess whether they are able to corrupt the transmission spectrum we infer for the planet.We used a flexible definition for spots as surface features that could be either cool (dark patches of magnetically-suppressed convection) or hot (bright patches like faculae or plage).Given all available data and a Poissonbased model for rotational variability, we estimated that the total spot covering fraction on GJ 1132

Figure 1 .
Figure1.A comparison of the average stellar spectrum extracted for each visit.We find no evidence of significant contamination by other sources as there is very little deviation in flux levels of the star across the five visits.The 1.15µm and 1.63µm are marked with dash lines for comparison with Figure2.

Figure 3 .
Figure 3. Broadband light curves for each WFC3 transit of GJ 1132b, with each visit plotted left to right, including the raw light curves with transit and ramp-like systematics (first row), the transits with the RECTE systematic model divided out (second row), the systematics with the BATMAN transit model divided out (third row), and the final residuals (fourth row).Outlying points are noted as light gray points.

Figure 4 .
Figure 4.Each of the 22 spectroscopic light curves plotted with the top to bottom being smallest to largest wavelength bins.Point colors correspond to one of the five visits.Light curves are shown with systematics removed and with the transit model created from the weighted-average transit depth of the five visits (left), and residuals from this shared model are shown (right).

Figure 5 .
Figure5.A transit depth comparison between the RECTE (green) and divide-white (purple) methods for each wavelength bin.Overall, transit depths vary less than 2σ for each bin (except for the second bin) with an average discrepancy of less than 1σ.

Figure 6 .
Figure 6.TESS folded light curve of GJ 1132b's transit.Top panel includes the detrended light curve points (gray), binned points for clarity (blue), and the best-fit model in light blue.Bottom panel shows the residuals from this fit.

Figure 7 .
Figure 7. Searching for TTVs in the mid-transit times of GJ 1132b yielded scatter around a flat line.Points are color-coded by telescope, and the HST uncertainties are inflated to match the 103-second exposure time of one observation point.

4.Figure 9 .
Figure 9.A comparison of the transit depths from each visit in their respective wavelength bins.The weighted-mean transit depths between the five visits are included in black.All bins demonstrated <2σ

4. 1 . 2 .Figure 10 .
Figure 10.A comparison between the average transit depths in Swain et al. (2021) (blue points), in Mugnai et al. (2021) (red points), and this work (black points).All utilize the same HST/WFC3 data, though Swain et al. (2021) determine very different transit depths, notably at shorter wavelengths.Transit depths derived in this work and in Mugnai et al. (2021) are statistically the same with similar uncertainties as well.
Southworth et al. (2017) observed 9 broadband transits with the GROND instrument on the 2.2m MPG telescope.They measure larger transit depths in both the z and K bandpasses, which suggest large H 2 O or CH 4 absorption in the atmosphere of GJ 1132b.Due to their large transit depth uncertainties and disagreement with other more precise data sets, we do not include theSouthworth et al. (2017) broadband points in any of our analyses.•Diamond-Lowe et al. (2018) used the LDSS3C spectrograph on the Magellan Clay telescope to observe the optical spectrum of GJ 1132b from 0.7 to 1.04 µm.They found a featureless spectrum and rule out a 1× Solar metallicity atmosphere with 3.6σ confidence.However, comparing their transit depths to the WFC3 spectrum depicts a slight offset, with their optical transit depths appearing smaller than those in the infrared.A close inspection shows that the broadband transit depth reported in Diamond-Lowe et al. ( Figure 11.A comparison between the WFC3 transit depths (black points) and the spectroscopic (dark blue) and broadband (light blue) transit depths from Diamond-Lowe et al. (2018).While the spectroscopic depths appear smaller than the WFC3 values, the broadband depth is statistically the same.

(
Kempton et al. 2017) for atmospheric compositions of 1×, 100×, 300×, 600×, and 1000× Solar metallicities by volume (Figure12) as well as atmospheres composed of 100% H 2 O, CO 2 , and O 2 (Figure13).These were all cloud-free and excluded any additional Rayleigh scattering beyond the gaseous species present in the atmosphere.Model pressure-temperature (P-T) profiles used the custom double-gray profiles (one opacity for short-wave radiation and another opacity for long-wave radiation) fromDiamond-Lowe et al. (2018).SeeDiamond-Lowe et al. (2018) for additional details on the calculations of these profiles.Stellar mass and radius were adopted fromBonfils et al. (2018).

Figure 14 .
Figure 14.The 0.7-4.5 µm transmission spectrum of GJ 1132b, combining all available archival data sets with as-published transit depths.The same atmospheric models from Figure 12 are included for comparison.

Figure 15 .
Figure 15.Similar to Figure 14, but with an offset of 0.000205 added to the Diamond-Lowe et al. (2018) spectroscopic transit depths.The addition of archival data strongly supports a featureless transmission spectrum from 0.7-4.5 µm for GJ 1132b.

DataFigure 17 .
Figure 17.MEarth out-of-transit light curve of GJ 1132, in the MEarth photometric bandpass (roughly 715-1000 nm), corrected for time-variable precipitable water vapor in Earth's atmosphere.Nightly bins are shown with a color intensity inversely proportional to the nightly variance, for two MEarth telescopes (tel13 in blue, tel16 in orange).The epochs of the five WFC3 visits are marked with vertical lines.
By using the gamma function Γ(x + 1) instead of x! in the denominator, we allow for non-integer values of the number of visible spots.Because we keep f 1 as a free parameter, this Poisson prior does not actually place new constraints on f and ∆f (t).Rather, it does allow us to connect them to the typical spot size f 1 that would best explain the amplitude of rotational variability coming from random Poisson draws for the number of starspots on the Earth-facing side of the star.The inferred typical spot size hinges on the assumption of how starspots are arrayed on the surface.Rackham   et al. ( ),Fang et al. (2016), andMorris et al. (2019) starspot samples.

Figure 18 .
Figure 18.A visualization of the posterior probability distribution for GJ 1132's starspot properties, including data used in the fit (left panels, with 20 random sample curves) and parameter samples (right panels, with 1000 random samples).Curves and points are colored from red to blue by T spot /T unspot , with paler colors closer to 1. Dashed lines highlight useful reference values in each parameter distribution.A fitted offset was applied to the LDSS3C depths, as explained in the text.The code used to generate this plot is available online (</>).
spectroscopic transit depths from Diamond-Lowe et al. (2018) demonstrate a significant offset from the other data sets due to possible instrumental effects.Fitting an offset to the Diamond-Lowe et al. (2018) points, we find that the majority of data sets (besides Southworth et al. 2017; Swain et al. 2021) rule out <300× Solar metallicity by volume atmospheric compositions with >5.6σ confidence.

Table 2 .
The R p /R s for each wavelength bin across each visit.We include the inverse-variance weighted average of the visits in the last column.To reproduce the transit depths plotted throughout this paper, we use (R p /R s ) 2 as the transit depth and 2(R p /R s )σ Rp/Rs to calculate the transit depth uncertainty.
data point.To more cautiously capture this uncertainty, we inflate the HST mid-transit uncertainties by a factor of 25 in order to match the 103 second exposure time of one point.Using the BJD TDB

Table 3 .
Individual HST and TESS Sector 9 Mid-Transit Times.
2.2 amu, 100% H 2 O is 18 amu, O 2 is 32 amu, and CO 2 is 44 amu.From these values, we calculate that 8.9 amu could represent an atmosphere mixed with 44% H 2 O and 56% H/He by volume.Model Figure 13.Same as Figure represents the absorption by the planet's atmosphere making the planet appear