Ground-based Optical Transmission Spectroscopy of the Small, Rocky Exoplanet GJ 1132b

Four years of transit data from the Kepler mission showed us that terrestrial planets are common around low-mass stars (Dressing & Charbonneau 2013, 2015; Gaidos et al. 2016). The Kepler data set also led to theories suggesting that some small planets retain hydrogen and helium envelopes from formation, comprising a small fraction of their total masses (Wolfgang & Lopez 2015). These H/He envelopes are subsequently sculpted by incident extreme ultra-violet (EUV) and X-ray radiation from the host stars, which, in the absence of a strong planetary magnetic field, drives atmospheric escape (Ehrenreich et al. 2015).

M dwarfs have extended pre-main-sequence phases (Baraffe et al. 2015) and remain chromospherically active on long timescales (Newton et al. 2017), so it is possible that terrestrial planets orbiting M dwarfs have been stripped of any primordial atmospheres early on (Lopez & Fortney 2013; Luger & Barnes 2015). For instance, the terrestrial planets TRAPPIST-1b and c orbiting an ultracool dwarf do not exhibit transmission spectra consistent with a cloud-free low mean molecular weight atmosphere at the level of ≥10σ confidence (de Wit et al. 2016). TRAPPIST-1d, e, and f also do not exhibit evidence for such atmospheres at the level of ≥4σ confidence (de Wit et al. 2018). We might expect a similar result for other small planets in close-orbits around cool stars.

In this work, we use ground-based observations to investigate the idea that terrestrial exoplanets orbiting M dwarfs do not possess low mean molecular weight atmospheres. We focus on the terrestrial exoplanet GJ 1132b (1.2 R_⊕, 1.6 M_⊕) orbiting a M4.5V dwarf that is 12 parsecs away from the solar system. The radius and mass of GJ 1132b are consistent with an iron and silicate composition similar to that of Earth and Venus (Berta-Thompson et al. 2015).

The surface gravity and estimated atmospheric temperature of GJ 1132b mean that a solar composition, hydrogen-dominated atmosphere might be detectable with ground-based instrumentation. Though we are looking for the signature of a low mean molecular weight atmosphere, hydrogen itself is not a strong absorber, making it difficult to detect via transmission spectroscopy. Instead, we assume the atmosphere to be well-mixed and search for tracer molecules like water (H₂O) or methane (CH₄), which have large absorption cross sections in the visible to near-infrared wavelengths.

Understanding the nature of terrestrial exoplanet atmospheres will bolster efforts to constrain planet formation and atmospheric evolution, and ultimately inform our search for biosignatures on other worlds. We do not expect life as we know it to exist on the highly irradiated surface of GJ 1132b, but understanding the atmospheres of hot, rocky planets will contextualize an eventual search for life on cooler, habitable zone exoplanets.

Though our current sample size of terrestrial exoplanets is small, it is important to understand them in the context of the well-studied solar system inner planets. Whether a terrestrial exoplanet resembles Earth or Venus or Mercury has vast implications for its formation history and life-hosting capabilities. Still more intriguing is the chance to uncover terrestrial planets with compositions and characteristics unseen in the solar system (e.g., Morley et al. 2017).

In Section 2, we describe our observations of GJ 1132b in transit. In Section 3, we describe our customized data-reduction pipeline, and in Section 4, we describe our customized data analysis pipeline. We present the results of this work in Section 5. We discuss the implications of ground-based investigations of terrestrial planet atmospheres in Section 6 and conclude with Section 7.

A joint program between Harvard and MIT (PIs Diamond-Lowe and Berta-Thompson, respectively) to observe transits of GJ 1132b received eight nights on the Magellan II (Clay) Telescope with the LDSS3C⁵ multi-object spectrograph at Las Campanas Observatory (Stevenson et al. 2016). Of the eight observing opportunities, we observed five transits of GJ 1132b and lost the remaining three nights to clouds and high winds. The details of our observing program are presented in Table 1.

GJ 1132 (V = 13.49, K = 8.322) is an M4.5V star (Berta-Thompson et al. 2015). In the 4' field of view of LDSS3C there are no stars of comparable magnitude or spectral type, so we opted to simultaneously observe nine comparison stars, which we later used to remove telluric effects from the GJ 1132 spectrum. Of these comparison stars, one was brighter than GJ 1132 but it saturated our detector and we were not able to use it in our analysis.

Our LDSS3C masks include slits for GJ 1132 and the nine comparison stars. At the time of our observations there was a background star 7.3'' away form GJ 1132; because GJ 1132 is a high proper motion star this separation will change over time and future observers of GJ 1132 should account for this. We oriented our mask such that the background star did not contaminate the dispersed spectrum of GJ 1132. We cut our slits 10'' in width to avoid slit losses and 20'' in length to provide sky background with which to perform our subtraction (Bean et al. 2010). We also cut identical masks with 1'' wide slits, which we used to take wavelength calibration arcs during the afternoon prior to each observation.

We set the detector binning to 2 × 2 and the readout speed to Fast (the LDSS3C user manual says this will give a 13 s readout time but we found it to be 16 s). We set the gain to Low, which, along with the readout speed, gives a gain of 0.6 ADU/electron. With our observation mask, we took biases, darks, quartz flat fields, and a mask image with which to align our stars in the slits during observations. With our 1'' mask, we took helium, neon, and argon arcs so that we could determine a wavelength solution for each dispersed stellar spectrum. Both during calibrations and observations we kept every detector pixel that we used to perform our analysis below 53000 ADU. As stated in the LDSS3C user manual and corroborated by the Las Campanas Observatory instrument specialists, the full pixel well is 65536 ADU, but past 53000 ADU, the detector stops counting photoelectrons linearly.

We chose to use the VPH-Red grism, which provides a wavelength coverage of 640–1040 nm with a central wavelength of 850 nm and a linear dispersion of 0.1175 nm/pixel (Stevenson et al. 2016). The VPH-Red grism has a higher resolution than the VPH-all grism, as well as a higher throughput at redder wavelengths. Using the VPH-Red grism allowed us to take longer exposures without saturating the detector, while also focusing on those wavelengths where GJ 1132 is brightest.

We took 13 s integrations and achieved a duty cycle of 45%. The VPH-Red grism introduces order contamination onto the detector, which we mitigated with the OG590 order-blocking filter as advised in the LDSS3C user manual. This filter blocks spectral contamination from higher spectral orders but produces internal reflections. (Stevenson et al. 2016 noted this contamination and decided against using the OG590 filter.) After inspecting the calibration arc frames during the day, we decided that the OG590 contamination was less problematic than the higher-order line contamination. We therefore used the OG590 filter during observation and also while taking our calibration images.

We note that our first night of observation (data set number 1 in Table 1) differed from the rest for two reasons. First, we neglected to use the OG590 order-blocking filter, which is why we exposed for 12 s on this night instead of 13. In spite of this, the order contamination was not drastic, as GJ 1132 emits few photons blueward of 700 nm. Second, we used a slightly different mask. The first amplifier (C1) of LDSS3C's CCD has several columns of bad pixels that over-lapped with one of our comparison stars. We cut a second, identical mask with the slits slightly shifted in order to avoid the bad pixels. We did not end up using this comparison star because the bad pixels near it saturated and leaked light into its dispersed spectrum. For consistency, we exclude this comparison star from all five data sets when performing the analysis.

For all five of our data sets, we acquired at least one transit-durations's worth of out-of-transit baseline both before and after the transit event with which to estimate the basline flux and correct for correlated noise in the data.

We transform our raw Magellan/LDSS3C images into one-dimensional (1D) stellar spectra by running them through our custom Python pipeline, mosasaurus.⁶ With this pipeline, we perform basic CCD processing on pairs of FITS images from the two amplifiers on LDSS3C. After subtracting 1D biases estimated from the amplifiers' overscan regions, we stitch images together into full frames, using the amplifiers' reported gains to convert from ADU to electrons. We create median-stacked two-dimensional bias and dark exposures that we subtract from all quartz flat and science exposures, to remove the baseline level of the readout electronics and the (very small) dark current accumulated during all exposures.

To identify and mitigate cosmic ray contamination, we compare each image to the 10 closest images in time. For each pixel, we calculate the median absolute deviation from the median (MAD), and flag any upward outliers that exceed 10× MAD as cosmic ray hits. We replace the flux value in the pixel affected by cosmic rays with its median value from the immediately surrounding exposures. We keep track of which pixels have been modified in this fashion, so they can be masked out of later analysis stages if so desired.

We cut out a 60 × 2048 pixel region around each of our dispersed spectra. We cut out corresponding regions using the same pixels on our quartz flat and arc images. The spectra recorded on the detector are curved slightly (by about 10 pixels over the entire thousand-pixel chip). We fit a quadratic function to the spectral trace, which maps where the peak flux is in each column in the cross-dispersion direction. We calculate the full-width-half-maximum (FWHM) in each column. To create an extraction window, we extend by three times the FWHM from the centroid in the cross-dispersion direction (Figure 1), allowing the extraction width to vary with wavelength. We create a range of extraction window sizes for each stellar spectrum, for later comparison. These apertures remain fixed with respect to the detector; they do not move to follow the slight motion of the spectral trace throughout the night (1–3 pixels).

Zoom In Zoom Out Reset image size

We use an interactive tool to plot the science images with an overlayed extraction window, to inspect the aligned extraction window containing the stellar flux, and to set custom sky subtraction regions uniquely for each star (Figure 1). For the GJ 1132 field, we determine at this stage that several comparison stars are unusable—the bright one that saturated the detector and four others that turned out to have multiple stars clustered together in the slit. Having multiple stars in a single slit is problematic as we would have to combine their spectra in a large extraction window, which leads to a poor estimate of the sky background and the blending of different wavelengths from different stars in the same extracted spectrum. We end up with four comparison stars for our analysis. Though we use the same comparison stars for each night of data, the extraction windows may vary from night to night. This is because the seeing conditions on a given night influence the PSF of the stars on the detector. We therefore stand to benefit from using different extraction windows for each star for each data set.

We median-pass filter quartz flat exposures, taken through the same wide slits as our science data, by dividing each pixel by the median of the 20 × 100 pixels surrounding it. We then divide each spectrum region in the time-series by this filtered quartz flat to correct for the intrinsic pixel-to-pixel inconsistencies of the detector. We create a 1D stellar spectrum from each flat-fielded stellar region by summing up the flux in each column in the cross-dispersion direction, accounting for partial pixels at the edges of the extraction aperture. We create a 1D sky-background spectrum from each flat-fielded stellar region by fitting a two-degree polynomial to each column in the cross-dispersion direction outside the extraction window and then summing over the column. We then subtract the sky-background spectrum from the stellar spectrum.

We tested an optimal extraction routine as outlined by Horne (1986). We find that this method makes at most a 10% improvement in signal-to-noise for the faintest comparison stars, but does not improve signal-to-noise for GJ 1132 or the brighter comparison stars. Because the fainter comparison stars have a proportionately small influence on the resulting light curve, we use the extraction method outlined above and not the optimal extraction routine.

We use the He, Ne, and Ar arcs taken during calibration to develop a rough wavelength solution for each star. The LDSS3C user manual provides a wavelength solution for the VPH-Red grism that gives the pixel position of prominent features in the He, Ne, and Ar spectra. Using a customized graphical user interface we match up the features in the provided wavelength solution to those in each arc, corresponding to our stellar spectrum regions, and create a polynomial wavelength solution for each star. In practice, this works better for some stars than others, but it generally lines up the spectra with each other to within 5 nm.

We then choose one exposure of one star as a basis against which to cross-correlate all of the exposures of all the stars in a given data set. We use five prominent features in the spectra in order to perform the cross-correlation: the O₂ doublet (760.5 nm), each line of the Ca triplet (849.8, 854.2, and 866.2 nm), and the forest of water lines (about 930–980 nm). We note that the Ca triplet is not a telluric feature and so may be risky to use when calibrating the spectra. In this case, all of the stars we observe are in the Sun's local moving group, and any Doppler shifting of the Ca lines are not detectable at the velocity dispersion of the LDSS3C spectrograph and the VPH-Red grism (about 165 km s⁻¹ pixel⁻¹). Given the small field-of-view of the instrument, we are not concerned about the different lines-of-sight to each star.

This process reveals that there is both a shifting and stretching of the spectra over the course of the observations. For instance, in data set number 1, the difference between the positions of the O₂ doublet and the water line forest in the GJ 1132 spectrum increases by a pixel from the start of the observation relative to the end. We use this information to apply a second wavelength solution for each spectrum in each exposure such that they are lined-up with one another in wavelength space to within 0.35 nm across all stars and the entire night. This is the final step in achieving 1D spectra that we can use to make our light curves.

To perform this analysis, we constructed the code detrendersaurus.⁷ Though it is not generalized for data sets other than LDSS3C multi-object spectroscopy, the code is fairly modular and some routines may be useful to others performing similar analyses.

4.1. Analyzing Transits Separately

GJ 1132 is brighter than the four comparison stars (Figure 2, Table 2). We therefore create our light curves by summing up the flux from the comparison stars and dividing the GJ 1132 spectrum by the summed comparison spectrum for each point in the light curve. GJ 1132 is still brighter than the summed flux of the four comparison stars, so we are limited by the photon noise of the summed comparison star flux.

Zoom In Zoom Out Reset image size

Table 2.  Stars Used in This Work

Star	R.A.	Decl.	Flux/
			GJ 1132 Flux
GJ 1132	10:14:50.09	−47:09:17.5	1.0
Comp A	10:14:57.51	−47:05:39.9	0.35
Comp B	10:14:58.22	−47:09:35.1	0.14
Comp C	10:15:05.74	−47:07:43.9	0.11
Comp D	10:15:16.26	−47:06:44.3	0.11

Note. The relative flux column indicates the full wavelength-band-integrated flux of each star relative to that of GJ 1132. The comparison star labels in this table correspond to those in Figure 2. GJ 1132 is a high proper motion star. Positions are given for an epoch of 2016.3.

Download table as:  ASCIITypeset image

We detrend our light curve using decorrelation parameters that either have the same values for all the stars (e.g., airmass) or are associated with GJ 1132 (e.g., width of the spectral trace). The parameters that are unique to each star have similar values for all stars in the data set but because we detect the most photons from GJ 1132 its decorrelation parameters have higher signal-to-noise ratios. We create white light curves for each data set and also bin the light curves from each data set into narrow wavelength bands for the purpose of atmospheric characterization. We restrict our analysis to the wavelength range common to GJ 1132 and the four comparison stars, which is 700–1040 nm.

We determine which linear combination of decorrelation parameters are necessary to remove the effects of correlated noise (discussed in greater detail in Section 4.2.1). In a given data set, we choose a single 20 nm wavelength bin without any prominent stellar, telluric, or atmospheric features (we use 830–850 nm) and calculate the Bayesian Information Criterion (BIC) value for every combination of possible decorrelation parameters. We check that there is no correlation with the sky-background, as this would imply that we are not properly removing the sky background during extraction. Once this process is done for all five data sets, we take the union of all the best decorrelation parameters and marginalize over them in all wavelength bins in all data sets. A list of these parameters, what vectors they depend on, and how they are derived can be found in Table 3.

Table 3.  Decorrelation Parameters Used to Model Data Systematics

Parameter	Vector	Description
Airmass	t	average airmass of the field
Rotator Angle	t	instrument rotator angle
Width	t, star	median width across wavelengths of the stellar trace in the cross-dispersion direction
Stretch	t, star	wavelength solution coefficient associated with spectrum stretching in the dispersion direction
Peak	t, λ, star	brightness of the brightest pixel in the cross-dispersion direction
Normalization	t	unit array

Note. All parameters are functions of time t. They can also vary by wavelength λ and by star. For all parameters that are star-dependent, we use the values associated with GJ 1132 as it has the highest signal-to-noise ratio.

Download table as:  ASCIITypeset image

From the results of a Levenberg–Marquardt minimizer, we run a makeshift Bayesian test in order to determine whether our five nights of data should be analyzed separately or taken together in a joint fit. For each 20 nm bin in each of the five data sets, we compare the resulting χ² value for a fit in which the transit depth is allowed to vary to one in which the transit depth is fixed to a inverse-variance weighted depth derived from the five nights. We account for the change in the number of fitted parameters between these two scenarios. We find that the χ² values for the case of the fixed transit depth in a given wavelength bin can be higher, lower, or identical to the case where the transit depth parameter is allowed to vary. In other words, fixing the transit depth does not provide a uniformly worse fit. We therefore decide to fit the five nights of data jointly, allowing the transit parameters to be shared across all nights.

4.2. Analyzing Transits Jointly

4.2.1. Levenberg–Marquardt Fits

In analyzing the transits jointly, we must account for the different uncertainties associated with the individual data sets, as well as clip outlying data points. We use a three-step Levenberg–Marquardt process to settle on initial guesses for our parameters to use in a dynamic nested sampler, which will be discussed in further detail in Section 4.2.2. To run our Levenberg–Marquardt fits, we employ the open-source lmfit package (Newville et al. 2016).

In the first pass at the Levenberg–Marquardt fit, we build a linear model unique to each night of data following the formula

$$\begin{eqnarray}&&{ \mathcal M }(t)={ \mathcal S }(t){ \mathcal T }(t)\end{eqnarray} \tag{ 1 }$$

where ${ \mathcal S }(t)$ is the systematics model and ${ \mathcal T }(t)$ is the transit model. The systematics model ${ \mathcal S }(t)$ can further be broken down to

$$\begin{eqnarray}&&{ \mathcal S }(t)=1+\displaystyle \sum _{n=1}^{N}{a}_{n}{p}_{n}(t)\end{eqnarray} \tag{ 2 }$$

where N is the number of decorrelation parameters used in the fit, a_n are the coefficients we are fitting for, and p_n are the arrays of decorrelation parameters that describe the correlated systematics in the data, which are all functions of time. For decorrelation parameters that are functions of wavelength, we sum over wavelength space corresponding to the wavelength bin we are working in. The decorrelation parameters are either common to all stars (airmass and rotator angle) or are taken from the GJ 1132 spectral extraction (width, stretch, peak).

We build the transit model ${ \mathcal T }(t)$ using the open-source batman code (Kreidberg 2015) and feed in the free transit parameters. The transit parameters that can be shared across the five data sets are the planet-to-star radius ratio ${R}_{p}/{R}_{* }$, period P, inclination i, scaled orbital distance $a/{R}_{* }$, and uncorrelated quadratic limb-darkening coefficients $2{u}_{0}+{u}_{1}$ and ${u}_{0}-2{u}_{1}$ as used by Holman et al. (2006). The residuals that we calculate from dividing our light curves by the linear models are weighted by the calculated photon noise of each data set.

At this stage, we fix the uncorrelated quadratic limb-darkening coefficients to values derived from the Limb Darkening Tool Kit (ldtk), an open-source package that takes in stellar parameters and uncertainties and calculates the limb-darkening coefficients in a given wavelength range based on the PHOENIX library of stellar models (Husser et al. 2013; Parviainen & Aigrain 2015). During the next stage of analysis (Section 4.2.2), we instead allow the uncorrelated quadratic limb-darkening parameters to vary within a prior.

In the second Levenberg–Marquardt fit, we calculate the MAD of the residuals and clip the 29 data points (for the white light curve) or ≤27 points (for the wavelength-binned light curves) that deviate by 5× the MAD. In the third Levenberg–Marquardt fit, we change the weighting from the calculated photon noise to the uncertainties we derive from each night's data as a result of our second fit. Levenberg–Marquardt fits with lmfit are inexpensive and quick but running a dynamic nested sampler can be expensive if the priors are too wide. Because we derive our sampling priors from the covariance matrix output by the Levenberg–Marquardt fit, we find it expedient to constrain the fit parameters as much as possible at this stage.

4.2.2. Dynamic Nested Sampling

Our joint fit comprises a minimum of 30 free parameters—the same six decorrelation parameters (Table 3) to fit for each of the five data sets. In addition to this, there can be free transit model parameters, like the transit midpoint for each night or the transit depth, which is shared between the five nights. Which transit parameters are free depends on whether or not we are performing a white light curve fit or a wavelength-dependent light curve fit. With so many free parameters traditional Markov Chain Monte Carlo (MCMC) ensemble samplers such as emcee (Foreman-Mackey et al. 2013) are slow and inefficient at exploring the parameter space (Huijser et al. 2015). We instead use the open-source dynamic nested sampling package dynesty⁸ (J. Speagle 2018, private communication) to estimate our posteriors.

The dynesty code samples each free parameter from 0 to 1 and so requires a prior transform function to map the outputs from the sampling onto the parameter space we want to explore. For all but the uncorrelated quadratic limb-darkening coefficients, we set uniform priors on the parameters used to model the systematic and transit portions of our models. When possible, we assume the same uniform priors for the transit model parameters as used by Dittmann et al. (2017a). Otherwise, we set uniform priors by taking the estimated 1σ uncertainties from the covariance matrix of our Levenberg–Marquardt fit and multiplying by 25 such that the prior bounds for each parameter are 25σ from the estimated parameter value. These wide uniform priors allow for an uninformed, broad parameter space for the sampler to explore.

Following the work of Berta et al. (2012), we place Gaussian priors on the uncorrelated quadratic limb-darkening coefficients $2{u}_{0}+{u}_{1}$ and ${u}_{0}-2{u}_{1}$. To determine what these Gaussian priors should be, we first get the quadratic limb-darkening coefficients in each wavelength bin from ldtk (Parviainen & Aigrain 2015). ldtk has an option to run an MCMC with the input stellar parameters and uncertainties in order to derive the limb-darkening coefficients. We use the samples from the MCMC to calculate arrays of uncorrelated parameters using the formulation $2{u}_{0}+{u}_{1}$ and ${u}_{0}-2{u}_{1}$, where u₀ and u₁ are the quadratic coefficients derived with ldtk. We calculate the median and standard deviation of these uncorrelated arrays and use these values to set the Gaussian priors. These Gaussian priors leverage our knowledge of stellar astrophysics without having to place complete faith in the accuracy of the stellar models.

Also following Berta et al. (2012), we achieve a χ² value of unity by including a rescaling parameter s (Equations (2) and (3) of that paper). We automatically marginalize over this during our dynamic nested sampling by modifying our log-likelihood function such that s is a multiplier of the theoretical uncertainty associated with each data point, including all terms that depend on s. Each data set has its own value of s associated with it. An s value of unity implies that we are reaching the theoretical photon noise limit with our fits, while a value less than unity implies an over-fitting of the model to the data.

4.2.3. White Light Curve

We jointly fit the white light curve of our five data sets and allow the time of mid-transit δ t₀ to vary for each data set, along with the shared parameters of the radius ratio ${R}_{p}/{R}_{* }$, period P, inclination i, scaled orbital distance $a/{R}_{* }$, and uncorrelated quadratic limb-darkening coefficients $2{u}_{0}+{u}_{1}$ and ${u}_{0}-2{u}_{1}$. In doing so, we leverage the five nights of data, which have the same transit model parameter values, except for the mid-transit time.

Where appropriate, we adopt the same priors as those quoted by Dittmann et al. (2017a, see our Table 4). The photometric bandpass of MEarth is not identical to the wavelength coverage of our white light curves, and so we use stellar models to set Gaussian priors on the uncorrelated quadratic limb-darkening coefficients, as described in Section 4.2.2.

Table 4.  White Light Curve Transit Model Parameter Priors

Parameter	Value	Priors
$\delta {t}_{\mathrm{0,1}}$ [days]	−0.0015	[−0.0075, 0.0044]a
$\delta {t}_{\mathrm{0,2}}$ [days]	−0.0017	[−0.0057, 0.0023]a
δt0,3 [days]	−0.0019	[−0.0108, 0.0069]a
δt0,4 [days]	−0.0016	[−0.0093, 0.0060]a
δt0,5 [days]	−0.0017	[−0.0060, 0.0027]a
${R}_{p}/{R}_{* }$	0.0493	[0.0081, 0.0904]a
P [days]	1.628925	[1.628744, 1.629116]b
i	88.68	[85, 90]b
$a/{R}_{* }$	16.54	[12, 20]b
$2{u}_{0}+{u}_{1}$	⋯	0.8756 ± 0.0128c
${u}_{0}-2{u}_{1}$	⋯	−0.3672 ± 0.0566c
s1,2,3,4,5	1	[0.01, 10]d

Notes.

^aUniform priors that are 25× the 1σ uncertainties taken from the lmfit covariance matrix, as described in Section 4.2.2. The δt₀ parameter is the offset from the calculated time of mid-transit (Equation (3)). ${R}_{p}/{R}_{* }$ is the planet-to-star radius ratio. ^bUniform priors taken from Dittmann et al. (2017a). P is the period, i is the inclination, and $a/{R}_{* }$ is the scaled orbital distance. ^cGaussian priors calculated with ldtk outputs, as described in Section 4.2.2, given as mean ± standard deviation. $2{u}_{0}+{u}_{1}$ and ${u}_{0}-2{u}_{1}$ are the uncorrelated quadratic limb-darkening parameters. In the Levenberg–Marquardt fits, these parameters are fixed to the ldtk outputs, but when sampling the parameter space with dynesty we use the Gaussian priors; dynesty does not require starting values as inputs. ^dWide uniform priors set by hand. Each data set has a rescaling parameter s as described in Section 4.2.2

Download table as:  ASCIITypeset image

For the time of mid-transit, we fit for an offset δt₀ from the calculated mid-transit time using the ephemeris T₀ given by Dittmann et al. (2017a):

$$\begin{eqnarray}&&\delta {t}_{0}={t}_{0}-({T}_{0}+{nP})\end{eqnarray} \tag{ 3 }$$

where n is the number of elapsed transits since the ephemeris transit, P is the period, and t₀ is the time of mid-transit for the n^th transit. We fit for the offset $\delta {t}_{0}$ as opposed to the mid-transit time itself in order to keep the model coefficients within a few orders of magnitude of each other. This is optimal for Levenberg–Marquardt fitting with lmfit.

In our full band-integrated white light curve fit from 700 to 1040 nm, we see significant features in the residuals. After experimenting with decorrelation parameters and wavelength clipping, we conclude that the deep water absorption bands redward of 920 nm are leaving imprints on the white light curve, suggesting changes in precipitable water vapor in Earth's atmosphere during some of our observations. The white light curves presented here do not include these problematic bands and are instead integrated from 700 to 920 nm.

At this stage, we investigate any transit timing variations by comparing our five derived mid-transit times to those quoted in the discovery paper (Berta-Thompson et al. 2015) and subsequent work with MEarth and Spitzer (Dittmann et al. 2017a) (Figure 3). The mid-transit times from the Spitzer data set reported by Dittmann et al. (2017a) are the BJD_OBS values taken from the Spitzer header files. We correct these values to BJD_TBD, which accounts for leap seconds. We use the values of P = 1.6289246 days and T₀ = 2457184.55804 days (Dittmann et al. 2017a) to calculate all times of mid-transit.

Zoom In Zoom Out Reset image size

We note that our times of mid-transit are consistently two minutes earlier than the predicted time. We have simultaneous transit observations with the MEarth telescope array, which do not agree with our transit times. The timing offset in this analysis could be due to some unexplored systematic in the instrument or the data reduction. We check our header-time conversions to BJD_TDB multiple times following Eastman et al. (2010), making sure to account for the exposure and readout times. As a test, we perform a simple data reduction using only polynomials and the batman transit light curve package (i.e., without the detrendersaurus pipeline) and were unable to derive transit times in agreement with those of MEarth and Spitzer.

This discrepancy does not affect our results with respect to the atmospheric analysis, as we fix the times of transit to the best-fit values when performing our atmospheric analysis, and the time of mid-transit does not affect the transit depth at the time resolution of our data.

We compare our derived values of the planet-to-star radius ratio ${R}_{p}/{R}_{* }$, period P, inclination i, and scaled orbital distance $a/{R}_{* }$ to those reported by Dittmann et al. (2017a) and find that our results are in agreement (Table 5). We present the raw white light curves, jointly fit white light curve, time-binned white light curve, and time-binned white light curve residuals in Figure 4.

Zoom In Zoom Out Reset image size

Table 5.  White Light Curve Derived Transit Model Values, Compared to Dittmann et al. (2017a)

Parameter	Value (this work)	Value (D17a)
δt0,1 [days]	$-{0.0015}_{-0.00023}^{+0.00022}$	⋯
$\delta {t}_{\mathrm{0,2}}$ [days]	$-{0.0017}_{-0.00017}^{+0.00017}$	⋯
$\delta {t}_{\mathrm{0,3}}$ [days]	$-{0.0017}_{-0.00054}^{+0.00055}$	⋯
$\delta {t}_{\mathrm{0,4}}$ [days]	$-{0.0016}_{-0.00027}^{+0.00027}$	⋯
$\delta {t}_{\mathrm{0,5}}$ [days]	$-{0.0017}_{-0.00018}^{+0.00017}$	⋯
${R}_{p}/{R}_{* }$	${0.0488}_{-0.0009}^{+0.0012}$	${0.0455}_{-0.0006}^{+0.0006}$
P [days]	${1.62893}_{-0.00013}^{+0.00013}$	${1.6289246}_{-0.0000030}^{+0.0000024}$
i [deg]	${88.54}_{-0.90}^{+0.90}$	${88.68}_{-0.33}^{+0.40}$
$a/{R}_{* }$	${15.91}_{-1.761}^{+1.236}$	${16.54}_{-0.71}^{+0.63}$
$2{u}_{0}+{u}_{1}$	${0.876}_{-0.012}^{+0.012}$	⋯
${u}_{0}-2{u}_{1}$	$-{0.371}_{-0.056}^{+0.055}$	⋯
${s}_{1}$	${4.27}_{-0.18}^{+0.19}$	⋯
${s}_{2}$	${2.73}_{-0.12}^{+0.12}$	⋯
${s}_{3}$	${6.08}_{-0.26}^{+0.28}$	⋯
${s}_{4}$	${5.24}_{-0.22}^{+0.24}$	⋯
s5	${2.89}_{-0.12}^{+0.13}$	⋯

Download table as:  ASCIITypeset image

4.2.4. Wavelength-binned Light Curves

We investigate the atmosphere of GJ 1132b by creating a transmission spectrum. We divide our light curves into 20 nm wavelength bins and jointly fit for ${R}_{p}/{R}_{* }$ and the uncorrelated quadratic limb-darkening coefficients $2{u}_{0}+{u}_{1}$ and ${u}_{0}-2{u}_{1}$, along with the systematic parameters for each respective data set. We fix the times of mid-transit t₀ for each night to the values determined from the white light curve fit. We fix the values of P, i, and $a/{R}_{* }$ in our binned wavelength fits to those reported by Dittmann et al. (2017a) as these are derived from a higher resolution Spitzer time-series.

Our joint fit produces single values for the radius ratio ${R}_{p}/{R}_{* }$ and the uncorrelated quadratic limb-darkening coefficients $2{u}_{0}+{u}_{1}$ and ${u}_{0}-2{u}_{1}$ for each wavelength bin, but each of the five data sets has its own linear fit to the systematics in the light curve. In order to make more meaningful comparisons between the systematic parameters in a given data set, we scale each of them by subtracting off the mean value and then dividing by the standard deviation. This ensures that all of our systematic parameters are on the same relative scale and so comparing their fitted coefficients describes the relative importance of each parameter to the fit (Figure 5).

Zoom In Zoom Out Reset image size

In Figure 6, we present our light curves after dividing out the systematic models for each data set. The wavelength boundaries, rms, transit depth, and median factor of the expected photon noise limit for each wavelength bin are given in Table 6. According to Figure 2 of Stefansson et al. (2017), our observations of GJ 1132 are limited by the photon noise so we did not estimate the scintillation noise for the analysis. Including scintillation noise would not change the resulting transit depths but it would decrease our values in the final column of Table 6.

Zoom In Zoom Out Reset image size

Across the 17 wavelength bands we achieve a median transit depth error of 90 ppm. We compare this to 80 ppm for two GJ 1132b transits with the Spitzer 4.5 μm channel and 55 ppm with 25 MEarth transits in its photometric band (Dittmann et al. 2017a).

We present our transmission spectrum in Figure 7 and compare it to two sets of four model transmission spectra generated by the Exo-Transmit open-source code (Kempton et al. 2017). As inputs we use custom double-gray temperature–pressure profiles and associated equation-of-state files as well as GJ 1132b's surface gravity and radius at 1 bar of atmosphere and GJ 1132's stellar radius (Miller-Ricci et al. 2009; Miller-Ricci & Fortney 2010).

The 1 bar planet radius is smaller than the transit radius by an amount that depends on the atmospheric composition, temperature, and gravity. As these values are uncertain, we allow the 1 bar planet radius to float in order to achieve the best transmission model fits to our data. Changing the 1 bar planet radius alters the amplitude of the model features as well as the overall depth of the model.

The Spitzer data from Dittmann et al. (2017a) can resolve the ingress and egress of a transit of GJ 1132b, so we adopt the stellar mass and radius quoted in that paper in order to create the temperature–pressure profiles and model transmission spectra.

Zoom In Zoom Out Reset image size

One set of four model transmission spectra assumes solar elemental abundances (dominant in H₂ and He) with metallicities that are 1, 10, 100, and 1000× solar by volume. In these solar composition atmospheres, the dominant sources of opacity that contribute to the transmission features are CH₄ and H₂O, with modest contributions from NH₃, H₂S, and K. Higher metallicity atmospheres have higher opacities, which strengthen the model features, but also higher mean molecular weights, which dampen the model features. These competing effects are the reason why they highest amplitude features are associated with the 10× solar metallicity model.

The other set of four model transmission spectra assume H₂ and H₂O atmospheres where H₂O makes up 1, 10, 50, and 100% of the atmosphere by volume. The solar composition models account for collision-induced absorption but the H₂/H₂O do not. Given how flat the model transmission spectra are this should not impact the results. All models assume a clear atmosphere (i.e., no aerosols).

We also compare the GJ 1132b transmission spectrum to a flat line at the inverse-variance weighted-average transit depth and to a linear fit to the transit depths. The wavelength bin-averaged values for the Exo-Transmit models are weighted by the recorded counts of a GJ 1132 spectrum across the same wavelength range. By using an observed spectrum of GJ 1132, we account for the difference in relative brightness of GJ 1132 as a function of wavelength, as well as the telluric features imprinted on the spectrum. Because our wavelength bins are so narrow, this weighting is virtually indistinguishable from a simple mean across the model wavelength bins. We use the wavelength bin-averaged values of the model transmission spectra to calculate the χ² values associated with the model fits to the measured transit depths.

Our results disfavor a clear, 1× solar metallicity atmosphere at 3.09σ (99.80%) and a clear, 10× solar metallicity atmosphere at 3.7σ (99.98%) confidence. We disfavor a 10% H₂O, 90% H₂ atmosphere at 3.5σ (99.95%) confidence. Our measured transmission spectrum is consistent with a flat line and with metallicities in excess of ∼10× solar or water abundances greater than ∼10%, for aerosol-free atmospheres.

We compare our results to those of other groups (Figure 8). Our spectrophotometric transit depths are in agreement with photometric transit depths from the MEarth survey and the Spitzer 4.5 μm bandpass (Dittmann et al. 2017a) but not in agreement with the photometric transit depths from the GROND multi-band imager (Southworth et al. 2017).

Zoom In Zoom Out Reset image size

6.1. Ground-based Detection of Terrestrial Exoplanet Atmospheres

Our data-reduction process highlights the difficulties of trying to detect terrestrial exoplanet atmospheres from the ground. The signal we are looking for is small (a transit depth of 0.24% and an atmospheric variation of 0.02%) and we are not able to reach the photon noise limit (Table 6). One question is whether more data could disfavor higher mean molecular weight atmospheres, or if less data was required to reach the same conclusions.

To answer this question we select a test-case 20 nm wavelength bin, from 830 to 850 nm, and jointly fit for ${R}_{p}/{R}_{* }$ and the uncorrelated quadratic limb-darkening coefficients $2{u}_{0}+{u}_{1}$ and ${u}_{0}-2{u}_{1}$, as we did in our analysis, using 1, 2, 3, 4, and 5 data sets in each fit. We add the data sets in order of decreasing signal-to-noise: first data set 2, then 5, 1, 4, 3. We record the error on the transit depth after each data set is added to the analysis. We compare these to the errors in transit depth from the first data set, scaled by the inverse of the square-root of the number of data sets included.

As shown in Figure 9, we require all five transits of GJ 1132b to rule out clear, low mean molecular weight atmospheres at 3σ confidence. Theoretically, eight transits are needed to rule out the highest mean molecular weight atmospheres we tested (1000× solar metallicity and 100% H₂O), though this is a minimum estimate, as we do not achieve the photon noise limit and therefore our error bars do not decrease by the square-root of the number of data sets included in the analysis.

Zoom In Zoom Out Reset image size

In the coming era of extremely large ground-based telescopes (ELTs) detecting and characterizing terrestrial exoplanet atmospheres may be in reach. For example, the Giant Magellan Telescope (GMT) will have a diameter of 24.5 m, compared to the 6.5 m diameter of Magellan Clay. This means that the GMT will receive about (24.5/6.5)² = 14.2 times the number of photons per observation. The science and instrument requirements for the GMT-Consortium Large Earth Finder (G-CLEF), an optical-band echelle spectrograph with a multi-object spectrograph setting and the first-light GMT instrument, suggest that with a single transit observation of GJ 1132b, GMT/G-CLEF would be able to rule out the high mean molecular weight atmospheres we tested in this analysis (Szentgyorgyi et al. 2014).

The caveat for all ELT observations is that reaching the photon noise limit from the ground will still be difficult. This difficulty is derived from the need for comparison stars for ground-based transmission spectroscopy. For the nearby systems that we are interested in, there are few comparison stars of similar spectral type and magnitude available. Increasing the field of view of spectrographs would allow for more and possibly better comparison stars, but, for the brightest stars, expanding the field of view sufficiently to include ideal comparisons will likely exceed the design capabilities of spectrographs. For such bright targets, it may be worth investing in ground-based high-resolution spectrographs (R > 100,000), which can make atmospheric detections without simultaneous observations of comparison stars (Snellen et al. 2013).

6.2. Theoretical Atmosphere of GJ 1132b

6.3. Searching for More Terrestrial Exoplanets

We investigate whether or not the small, rocky terrestrial exoplanet GJ 1132b possesses a low mean molecular weight (μ ∼ 2) atmosphere using ground-based telescopes and instrumentation to construct a transmission spectrum. Our analysis disfavors a clear, 10× solar metallicity and a clear 10% H₂O at high confidence. GJ 1132b likely possesses a high mean molecular weight or depleted atmosphere.

While we search for new terrestrial exoplanets, we should also continue to learn more about the GJ 1132b atmosphere. Obtaining transits with HST/WFC3 will allow us to confirm the results from this work, especially considering that space-based telescopes do not have to contend with telluric water features. Morley et al. (2017) suggest that GJ 1132b is the most amenable planet of its kind, currently known, for observation in secondary eclipse with JWST. Small, rocky exoplanets like GJ 1132b challenge our limits of detection and characterization but also present the most exciting opportunities for comparative planetology with the solar system terrestrial exoplanets, including Earth.

This paper includes data gathered with the 6.5 m Magellan Telescopes located at Las Campanas Observatory, Chile. We thank the contributors to the LDSS3C project, the telescope operators and staff at Las Campanas Observatory, and the writers and contributors of the open-source software used in this work. We also thank Mercedes Lopez-Morales, Robin Wordsworth, Dimitar Sasselov, Laura Kreidberg, Kevin Stevenson, and Josh Speagle for helpful comments and conversations. H.D.-L. recognizes support from the National Science Foundation Graduate Research Fellowship Program (grant number DGE1144152). Z.K.B.-T. acknowledges support from the MIT Torres Fellowship for Exoplanetary Research. E.M.-R.K. was supported by the National Science Foundation under CAREER grant No. 1654295 and by the Research Corporation for Science Advancement through their Cottrell Scholar program. This publication was made possible through the support of a grant from the John Templeton Foundation. The opinions expressed here are those of the authors and do not necessarily reflect the views of the John Templeton Foundation.

Software: lmfit (Newville et al. 2016), dynesty (github.com/joshspeagle/dynesty), batman (Kreidberg 2015), Exo-Transmit (Kempton et al. 2017), mosasaurus (https://github.com/zkbt/mosasaurus), detrendersaurus (https://github.com/hdiamondlowe/detrendersaurus).