Detection of the Atmosphere of the 1.6 M_⊕ Exoplanet GJ 1132 b

M dwarfs are bounteous throughout our Galaxy and host large numbers of planets (Cassan et al. 2012; Dressing & Charbonneau 2015). The myriad of planets orbiting M dwarfs is dominated by small and low-mass rocky bodies with masses and radii comparable to Earth's (Morton & Swift 2014; Gaidos et al. 2016), and many occur in multiple systems (Muirhead et al. 2015). Planets around M dwarfs are of particular interest because the dimness of the host stars means that their habitable zones (e.g., Kopparapu et al. 2013) are located at short orbital periods, which are much more accessible to observational study (Gillon et al. 2016). However, planets within the habitable zone of M dwarfs face additional challenges such as tidal locking and large incident fluxes of high-energy photons (Lammer et al. 2003; Cunha et al. 2015; Shields et al. 2016).

From an analysis of the M dwarfs observed by the Kepler satellite, Gaidos et al. (2016) found that there were on average 2.2 planets per star and that the radius distribution peaked at 1.2 ${R}_{\oplus }$. These objects are expected to have tenuous atmospheres, with those of radius below approximately 1.5–1.6 ${R}_{\oplus }$ being almost entirely rocky (Weiss & Marcy 2014; Rogers 2015). This raises a problem for the overarching goal of constraining the chemical compositions and atmospheric characteristics of rocky planets, because the observational signatures of their atmospheres are below the levels detectable with current facilities.

GJ 1132 is a nearby very low mass star that was recently found to host a transiting and potentially rocky planet (Berta-Thompson et al. 2015, hereafter BT15) of mass 1.6 ${M}_{\oplus }$ and radius 1.2 ${R}_{\oplus }$. Its proximity to the Sun (12.04 ± 0.24 pc; Jao et al. 2005) means that it is comparatively bright and therefore well suited to analyses aimed at constraining the properties of both the star and the planet. Schaefer et al. (2016) presented simulations of the interior and atmosphere of GJ 1132 b, finding that the atmosphere is likely tenuous and dominated by O₂. While significantly hotter than Earth, the planet is one of the coolest transiting planets of known mass and is therefore of great interest for comparative planetology. BT15 found that GJ 1132 A is an old (5 Gyr or more) and slowly rotating (rotation period of 125 days) M4.5 V star, making it representative of a large population of planet host stars expected to be found in the (relatively) near future (Cloutier et al. 2017).

Dittmann et al. (2016, hereafter D16) presented extensive photometry of the GJ 1132 system, comprising 21 transits observed with the MEarth telescopes and a 100 hr light curve from the Spitzer satellite. D16 presented revised properties for the system and searched for transit timing variations or additional transits that would be caused by the presence of a third body in the system.

In this work we present extensive simultaneous optical and near-infrared photometry of nine transits of the planet GJ 1132 b in front of the star GJ 1132 A. We use these data to redetermine the physical properties of the system, improve the fidelity of its orbital ephemeris, and construct a transmission spectrum of the planet. We then interpret the transmission spectrum using suites of theoretical spectra from two model atmosphere codes. We clearly detect the planetary atmosphere, but the data in hand are not able to resolve ambiguities in the relative contributions of different molecules to the atmospheric opacity.

At 1.6 ${M}_{\oplus }$ GJ 1132 b is by a substantial factor the lowest-mass planet with an atmosphere that has been observationally detected. The two other low-mass planets with claimed detections are 55 Cnc e (Winn et al. 2011), with a mass of 8.08 ± 0.31 ${M}_{\oplus }$ (Demory et al. 2016), for which Tsiaras et al. (2016) found atmospheric features that could most easily be explained by HCN opacity, and GJ 3470 b (Bonfils et al. 2012), with a mass of 13.7 ± 1.6 ${M}_{\oplus }$ (Biddle et al. 2014), which shows an enhanced radius in the blue attributable to Rayleigh absorption (Nascimbeni et al. 2013; Biddle et al. 2014; Dragomir et al. 2015).

Three other low-mass planets have been subjected to transmission spectroscopy that has failed to reveal atmospheric signatures, most likely due to the presence of clouds or of a high-metallicity atmosphere with a large mean molecular mass. They are GJ 1214 b (Charbonneau et al. 2009), with a mass of 6.26 ± 0.91 ${M}_{\oplus }$ (Anglada-Escudé et al. 2013), for which no atmospheric features have been detected (Bean et al. 2010, 2011; de Mooij et al. 2012; Kreidberg et al. 2014b); HD 97658 b (Dragomir et al. 2013), with a mass of ${7.55}_{-0.79}^{+0.83}$ ${M}_{\oplus }$ (Van Grootel et al. 2014), for which there was also a nondetection of the atmosphere (Knutson et al. 2014b), and GJ 436 b (Gillon et al. 2007), with a mass of 25.4 ± 2.1 ${M}_{\oplus }$ (Lanotte et al. 2014), for which Knutson et al. (2014a) found no atmospheric features.

Extensive observations of GJ 1132 were obtained in service mode using the GROND multiband imager (Greiner et al. 2008) mounted on the MPG 2.2 m telescope at ESO La Silla, Chile. This instrument acquires images simultaneously in four optical and three near-infrared passbands. A total of nine transits were observed, using the telescope-defocusing approach (Alonso et al. 2008; Southworth et al. 2009), with point-spread functions (PSFs) of typically 30 pixels in radius. A 10th transit observation suffered from a large systematic effect during transit, due to either weather or instrumental effects, so it was not included in our analysis.

The optical bands have response functions similar to those of the SDSS g, r, i and z filters (Fukugita et al. 1996). They are equipped with 2k × 2k CCDs that see (approximately) the same 54 × 54 field, with a plate scale of 0158 pixel⁻¹, and are each read out through two amplifiers. The near-infrared bands have 1k × 1k Rockwell HAWAII-1 detectors with a field of view of 10' × 10', which fully encompasses the optical field, at a plate scale of 06 pixel⁻¹.

Several phenomena affected the quality of our optical observations. First, GJ 1132 is observationally difficult for the GROND imager, due to its low ${T}_{\mathrm{eff}}$ and concomitant large variations in flux level through the optical wavelength region. The requirement for a common exposure time and focus level meant that a compromise had to be made between obtaining adequate count rates in the g band and avoiding overexposure in the z band. Three of the z-band observing sequences suffered from saturation effects that precluded the extraction of reliable photometry from these images. Second, an issue was encountered in the r and z bands attributable to electronic noise in the CCD controllers. This caused approximately 0.3% of the pixels on each image to be assigned count rates near the bias level, with the identities of the affected pixels varying at random with each image. When these coincided with the PSFs of the target or comparison stars, they caused a nonastrophysical drop in the number of counts detected from the star, increasing the scatter on the photometry. A nonrepeating effect such as this cannot be calibrated out, but a partial mitigation of the effect was achieved by detecting each pixel with anomalously low count rates and replacing its value with the mean level from the adjacent 8 pixels. Third, the small field of view of GROND meant that the r, i, and z bands suffered from a lack of decent comparison stars. Photometry was therefore performed against comparison stars that were significantly fainter than GJ 1132 itself. A brief observing log of the nine included transits is given in Table 1.

The optical data were reduced using the defot pipeline (Southworth et al. 2009, 2014), which performs aperture photometry using the aper routine from daophot (Stetson 1987) contained in the NASA astrolib library.⁷ The aperture positions and sizes were specified manually in order to yield data with the lowest statistical and systematic errors. The science images were not calibrated using bias or flat-field frames, as these tended to have little effect on the final light curves beyond a slight increase in the scatter of the data points. Differential magnitudes for each light curve were formed versus an optimal ensemble of comparison stars, while simultaneously fitting for and removing a quadratic trend of magnitude with time caused primarily by airmass variations. The time stamps were converted into the BJD(TDB) system using routines from Eastman et al. (2010). As a final step, we rescaled the error bars for each light curve to obtain a reduced ${\chi }^{2}$ of ${\chi }_{\nu }^{\,2}=1.0$ versus a fitted model (see below). The light curves are shown in Figure 1 and will be made available at the CDS.⁸

Zoom In Zoom Out Reset image size

The infrared data were reduced using standard IDL routines following the methods outlined by Chen et al. (2014). In brief, we constructed a master dark frame and a master flat field by median-combining the individual dark and sky flat-field images obtained before the science observations. From each science image we subtracted the master dark frame and divided by the normalized master flat field. No correction for nonlinearity was applied to the data because the counts were below the level at which linearity becomes important for the GROND infrared detectors. We obtained the light curves using aperture photometry routines to extract the flux of the target and several comparison stars. We also tried to correct for the odd–even readout pattern present along the x-axis, but found no improvement in the light curve, so we decided not to include the correction for this effect when obtaining our final light curves. We were able to obtain useful results for five transits in each band (Table 1); the remaining data sets suffered from correlated noise features that were larger than the transit depth.

The optical GROND data were combined into four sets, one for each of the g, r, i, and z passbands, and then each was modeled using the jktebop code (Southworth 2013). We fitted for the orbital inclination (i), time of midtransit (T₀), and the sum and ratio of the fractional radii (${r}_{{\rm{A}}}+{r}_{{\rm{b}}}$ and $k={r}_{{\rm{b}}}/{r}_{{\rm{A}}}$), where the fractional radii are those of the star and planet in units of the orbital semimajor axis (${r}_{{\rm{A}},{\rm{b}}}=\tfrac{{R}_{{\rm{A}},{\rm{b}}}}{a}$). The orbital period was fixed to the known value (see Section 4), and the orbit was assumed to be circular. A quadratic function versus time was applied to the magnitude values for each observed transit in order to propagate the uncertainty in this from the light-curve generation process.

Limb darkening was incorporated using each of five laws (see Southworth 2008) and with the nonlinear coefficient fixed. Fits were obtained with the linear coefficient either fixed or fitted in order to see how this changed the results. We adopted limb-darkening coefficients calculated using the phoenix model atmospheres by Claret (2004), and linearly interpolated the coefficients to the stellar temperature of ${T}_{\mathrm{eff}}=3270\pm 140$ K (BT15).

An additional complication arises owing to the presence of the faint nearby star USNO B1.0 0428–0265237 (Monet et al. 2003), which was separated from GJ 1132 by approximately 65 at the times of our observations. A small fraction of the flux from this star leaked into the aperture used to measure the counts of GJ 1132. From a simple PSF model we measured the amount of contamination as 1.2% ± 0.3% in the g band, 0.5% ± 0.2% in r, 0.3% ± 0.1% in i, and 0.2% ± 0.1% in z, where the error bars are very conservative. This effect was included as "third light" in the jktebop model following the approach in Southworth (2010). The best fits are plotted in Figure 2 in the case of fixed limb-darkening coefficients.

Zoom In Zoom Out Reset image size

We also modeled the two best light curves of GJ 1132 from BT15, which are the g- and i-band data from the PISCO imager. Their appearance in Figure 2 differs from that in BT15, leading us to investigate the discrepancy. We checked for correlation with the instrumental parameters supplied with the flux measurements (airmass, x-position, y-position, PSF width, and sky background level), finding linear Pearson correlation coefficients less than 0.13 in all cases. We therefore conclude that the difference in appearance is only because BT15 binned their data into 1.5-minute intervals before plotting it.

Shortly before our own manuscript was submitted, D16 presented an analysis of GJ 1132 based on observations of 21 transits with the MEarth telescopes and 100 hr with the Spitzer satellite using the 4.5 μm channel of the IRAC imager. While the original data are not available to us, we have been able to obtain binned light curves from both facilities by digitizing Figure 2 in D16. We used the scatter of the data around the best-fit model to assess the photometric precisions of the data sets. For the Spitzer data we used limb-darkening coefficients from Claret et al. (2012), and for the MEarth data we adopted the mean of the coefficients for the i and z bands used above. For the Spitzer light curve we found that fitting for any limb-darkening coefficients returned a net limb brightening, due to the scatter of the data during transit. The MEarth light curve was able to support the fitting of one limb-darkening coefficient. Results for both data sets are included in Table 2.

Table 2.  Parameters of the Fit to the GROND and PISCO Light Curves of GJ 1132 from the jktebop Analysis

Source	${r}_{{\rm{A}}}+{r}_{{\rm{b}}}$	k	i (deg)	${r}_{{\rm{A}}}$	${r}_{{\rm{b}}}$
GROND g band	0.102 ± 0.028	0.0535 ± 0.0060	85.4 ± 2.0	0.097 ± 0.026	0.0052 ± 0.0019
GROND r band	0.066 ± 0.019	0.0517 ± 0.0040	88.5 ± 2.1	0.063 ± 0.018	0.0032 ± 0.0012
GROND i band	0.070 ± 0.014	0.0514 ± 0.0023	87.9 ± 1.6	0.067 ± 0.013	0.0034 ± 0.0008
GROND z band	0.106 ± 0.028	0.0611 ± 0.0061	85.3 ± 1.9	0.100 ± 0.026	0.0061 ± 0.0021
PISCO g-band	0.117 ± 0.039	0.0633 ± 0.0086	84.5 ± 2.6	0.110 ± 0.036	0.0070 ± 0.0029
PISCO i band	0.086 ± 0.030	0.0541 ± 0.0050	86.6 ± 2.4	0.082 ± 0.028	0.0044 ± 0.0019
Spitzer  4.5 μm	0.090 ± 0.014	0.0496 ± 0.0010	86.3 ± 1.0	0.086 ± 0.013	0.0043 ± 0.0007
MEarth	0.064 ± 0.024	0.0492 ± 0.0026	88.5 ± 1.8	0.061 ± 0.024	0.0030 ± 0.0013
Final results	0.0814 ± 0.0072	0.05041 ± 0.00086	86.58 ± 0.63	0.0775 ± 0.0068	0.00397 ± 0.00042

Download table as:  ASCIITypeset image

Table 2 gives the best-fitting photometric parameters for each light curve, together with the weighted mean value for each parameter. Error estimates for the photometric parameters were obtained using Monte Carlo (Southworth et al. 2004) and residual-permutation simulations (Southworth 2008), and the larger value retained in each case. These were then inflated to account for any variation between the solutions adopting the five different limb darkening laws by adding in quadrature the largest difference in values of each parameter from solutions with the various limb darkening laws. Finally, the weighted mean values agree well for all parameters, indicating that the error bars are robust. In all cases we have chosen to adopt the solutions with all limb darkening coefficients fixed. We find that the alternative solutions with one fitted coefficient agree to well within the uncertainties, with the exception of the Spitzer data as discussed above.

There is a significant disagreement between the photometric parameters determined by us and those determined by BT15 and D16. An inspection of Table 2 shows that there is a correlation between the parameters ${r}_{{\rm{A}}}$ and i, whereby higher-inclination solutions give a smaller fractional radius for the star. BT15 and D16 obtain solutions with higher i than our own, resulting in lower ${r}_{{\rm{A}}}$ values of 0.0625 ± 0.0043 and ${0.0605}_{-0.0022}^{+0.0027}$, respectively. This can be explained in the case of BT15 because ${r}_{{\rm{A}}}$ is a proxy for the density of the star (Sozzetti et al. 2007; Southworth 2017), so this quantity was effectively fixed in the light-curve solutions by the imposition of external constraints on both the mass and radius of the host star.

Such a constraint was not imposed in the current work or D16, so the discrepancy between these two works remains. We have checked and ruled out the possibility that it could be caused by the treatment of limb darkening or third light, by quantifying the effect of varying our own treatment of these phenomena within reasonable limits. We also tried numerically integrating the fitted model (Southworth 2011) to match the cadence of the Spitzer and MEarth data, but found that this caused little change in the solutions. We also tried various initial values for the fitted parameters, including the values found by D16, but the fits to data converged on the same solutions as given in Table 2. Solution determinacy is therefore not an issue, and the cause of the discrepancy remains unidentified. It is most likely due to differences in the overall analysis procedure.

Many of our light curves have a relatively high scatter compared to the transit depth. We therefore combined all four light curves of each individual transit into a single data set and fitted them with jktebop while fixing the sum of the radii (${r}_{{\rm{A}}}+{r}_{{\rm{b}}}$) to the best fit from Table 2. This yielded nine measured times of midtransit, which are given in Table 3.

We augmented these timing measurements with the reference time of midtransit given by BT15 from a global analysis of their photometric data (2,457,184.55786 ± 0.00032) and with the 22 timings given by D16 from orbital cycle 91 onward. The timings were then fit with a straight line to yield the linear ephemeris:

$$\begin{eqnarray*}&&{T}_{0}=\mathrm{BJD}(\mathrm{TDB})\,2457184.55759(30)\,+\,1.6289287(18)\times E,\end{eqnarray*}$$

where the fit has ${\chi }_{\nu }^{\,2}=1.13$ and the uncertainties (given in brackets and relative to the preceding digit) have been increased by $\sqrt{1.13}$ to account for this. The scatter around the best fit gives no indication of deviations from a linear ephemeris, so we do not find any evidence for transit timing variations.

The physical properties of the GJ 1132 system were established by BT15 using the following steps. First, the measured trigonometric parallax (BT15; Jao et al. 2005) and Two Micron All Sky Survey JHK magnitudes were used to determine the K-band absolute magnitude of the star, from which its mass was found via the calibration by Delfosse et al. (2000). Second, an empirical calibration of mass versus density (Hartman et al. 2015) was used to find the density and thus radius of the star. Third, the transit light curves were fitted with the stellar density constrained to the value found in the second step, which is equivalent to constraining ${r}_{{\rm{A}}}$.

We are in a position to modify this approach to rely less on general empirical calibrations and more on data obtained for GJ 1132 itself. We adopted the stellar mass of 0.181 ± 0.019 ${M}_{\odot }$ from BT15 and used this, along with the photometric parameters measured in Section 3, to determine the physical properties of the system using standard equations (e.g., Hilditch 2001). This process in effect used ${r}_{{\rm{A}}}$ as a proxy for the density of the star (Seager & Mallén-Ornelas 2003), which, combined with its mass, yields its radius. The properties of the planet could then be determined relative to those of the star using k and the velocity amplitude of the star measured by BT15, 2.76 ± 0.92 m s⁻¹. The uncertainties were propagated by a Monte Carlo approach.

The measured physical properties of the GJ 1132 system are summarized in Table 4. It can be seen that we find a rather lower density for the star—measured directly from the light curves, whereas BT15 obtained their value from a calibration of stellar properties—which results in increased radii for both star and planet. The measured density of GJ 1132 b has decreased by a factor of two (1.8σ), which has a significant impact on the likely properties of this body.

Table 4.  Derived Physical Properties of GJ 1132 from the Current Work, in Comparison with Values Found by BT15 and D16

Quantity	Symbol	Unit	This work	BT15	D16
Stellar mass	${M}_{{\rm{A}}}$	${M}_{\odot }$	0.181 ± 0.019	0.181 ± 0.019	⋯
Stellar radius	${R}_{{\rm{A}}}$	${R}_{\odot }$	0.255 ± 0.023	0.207 ± 0.016	${0.2105}_{-0.0085}^{+0.0102}$
Stellar surface gravity	$\mathrm{log}{g}_{{\rm{A}}}$	cgs	4.881 ± 0.074	⋯	⋯
Stellar density	${\rho }_{{\rm{A}}}$	${\rho }_{\odot }$	${10.9}_{-2.4}^{+3.4}$	21.0 ± 4.3	${19.4}_{-2.5}^{+2.6}$
Planet mass	${M}_{{\rm{b}}}$	${M}_{\oplus }$	1.63 ± 0.54	1.62 ± 0.55	⋯
Planet radius	${R}_{{\rm{b}}}$	${R}_{\oplus }$	1.43 ± 0.16a	1.16 ± 0.11	1.130 ± 0.056
Planet surface gravity	${g}_{{\rm{b}}}$	m s−2	${7.8}_{-2.8}^{+3.4}$	11.7 ± 4.3	⋯
Planet density	${\rho }_{{\rm{b}}}$	g cm−3	${3.1}_{-1.2}^{+1.7}$	6.0 ± 2.5	6.2 ± 2.0
Equilibrium temperature	${T}_{\mathrm{eq}}\,^{\prime}$	K	644 ± 38	579 ± 15	⋯
Orbital semimajor axis	a	au	0.01533 ± 0.00053	⋯	⋯

Note.

^aThis value of R_b is averaged over all observed passbands. However, while interpreting the spectrum, we consider the z-band radius at 1.57 ± 0.05 ${R}_{\oplus }$ to be contributed by the planetary atmosphere, whereas the remaining data points are consistent with a continuum, the "surface" of the planet, at a bulk radius of 1.375 ${R}_{\oplus }$.

Download table as:  ASCIITypeset image

The main aim of the current work is to explore whether constraints on the atmospheric composition of GJ 1132 b can be obtained using multiband photometry. We therefore sought to determine the radius of the planet in each of the passbands for which we possess a light curve. It is important to remove sources of uncertainty common to all passbands, so that the significance of any relative variations between passbands can be assessed.

Following the approach of Southworth et al. (2012), we modeled the nine available light curves (GROND griz, GROND JHK, and PISCO gi) with ${r}_{{\rm{A}}}$, i, and the orbital period fixed to the final values determined above. We included as fitted parameters ${r}_{{\rm{b}}}$, the reference time of midtransit, and the quadratic function of magnitude versus time for each transit. We adopted quadratic limb darkening with the linear coefficient fitted and the nonlinear coefficient fixed to the values used above. As the ${r}_{{\rm{A}}}$ and i parameters are fixed, it is not only more tractable but also more important to fit for limb darkening in order to propagate its uncertainty and avoid biases arising from the use of theoretical values. Third light was included for the optical bands and constrained as in Section 3. For the JHK bands we neglected third light, as it is negligible at infrared wavelengths, but iteratively clipped points lying more than 3σ from the best fit in order to remove scattered data that was biasing the fitting process. The near-IR light curves and best fits are shown in Figure 3. We did not include the data from D16 in this analysis because imperfections in the digitization process could significantly affect our results. In addition, the wavelength resolution of the MEarth data is poor, and the Spitzer data cannot support the fitting of limb darkening.

Zoom In Zoom Out Reset image size

For each light curve we determined the best-fitting ${r}_{{\rm{b}}}$ and obtained its uncertainty via 1000 Monte Carlo simulations. The results of this process are given in Table 5, along with details of the passbands used, and shown in Figure 4. For the purposes of visualizing our results, we assumed that the PISCO passbands correspond to those from Fukugita et al. (1996). It can be seen that the g, r, and i passbands agree well, except for a discrepancy between our own g-band results and those from the PISCO g-band data from BT15. We suggest that our ${r}_{{\rm{b}}}$ value is more reliable as it is based on observations of nine transits, whereas the PISCO data cover only one transit. It is also apparent that the radius of the planet in z is significantly larger than that in the other passbands, the discrepancy being a maximum of 4.1σ versus our i-band value of ${r}_{{\rm{b}}}$.

Zoom In Zoom Out Reset image size

Table 5.  Values of ${r}_{{\rm{b}}}$ for Each of the Six Light Curves

Passband	Central	Band Full	${r}_{{\rm{b}}}$
	Wavelength (nm)	Width (nm)
GROND g	477	138	0.00382 ± 0.00011
GROND r	623	138	0.00402 ± 0.00009
GROND i	763	154	0.00386 ± 0.00006
GROND z	913	137	0.00446 ± 0.00015
PISCO g	⋯	⋯	0.00438 ± 0.00010
PISCO i	⋯	⋯	0.00396 ± 0.00007
GROND J	1230	410	0.00354 ± 0.00045
GROND H	1645	420	0.00324 ± 0.00044
GROND K	2165	570	0.00473 ± 0.00058

Note. The error bars in this table exclude all common sources of uncertainty and so should only be used to interpret relative differences in ${r}_{{\rm{b}}}$. To convert ${r}_{{\rm{b}}}$ to ${R}_{{\rm{b}}}$, we multiplied by 23,455.0, the orbital semimajor axis in units of ${R}_{\oplus }$. The central wavelengths and FWHM transmission are given for the GROND filters.

Download table as:  ASCIITypeset image

6.1. Testing the Robustness of the Results

The optical transmission spectrum is the most important result of the current work, so we have deployed a suite of tests to probe the reliability of the radius values and error bars found for each passband. These have been applied to the optical bands only, as it is here that the measurements are most precise and extensive, and because the optical bands are the most affected by issues such as limb darkening and contaminating light.

6.1.1. Sensitivity to Individual Transit Observations

A cursory inspection of Figure 1 reveals that the transit depth in each passband may vary between different transits. The most obvious indicators are an unusually shallow transit in g on 2016 March 26 and an apparently very deep transit in z on 2016 May 17. This possibility can be probed by determining a set of solutions in each passband, each one based on the full data set but with one transit omitted. We performed this analysis for the GROND data and plot the results in Figure 5. The PISCO data were not considered: they are unsuitable for this analysis because they cover only one transit.

Zoom In Zoom Out Reset image size

We found that the planetary radius measurements were not significantly affected by the omission of individual transit data sets: all results exist within the 1σ error bars found for our default solution. We conclude that our results are robust against the omission of individual light curves. This is a particularly powerful test in the current case, because of the large number of data sets obtained in each band, and makes us confident in the reliability of both the measured values and error bars of the optical transmission spectrum.

We also modeled every optical transit individually and show the resulting ${R}_{{\rm{b}}}$ values in Figure 6. The z-band values of ${R}_{{\rm{b}}}$ show a clear increase relative to those for other bands, peaking at a 4.4σ difference between i and z. Even if we reject the highest radius measurement, the last one in the z band, we obtain a difference in ${R}_{{\rm{b}}}$ between the i and z bands significant at the 3.1σ level. The measurements exhibited in Figure 6 are tabulated in Table 6.

Zoom In Zoom Out Reset image size

Table 6.  Measured Planetary Radius for Each Individual Transit in the GROND griz Bands, with Relative Error Bars, as Plotted in Figure 6

Observing Night	Measured Planetary Radius (${R}_{\oplus }$)
	g	r	i	z
2016 Jan 14	1.209 ± 0.154	1.393 ± 0.125	1.420 ± 0.062	1.459 ± 0.105
2016 Feb 14	1.475 ± 0.091	1.302 ± 0.110	1.422 ± 0.051	⋯
2016 Feb 24	1.567 ± 0.151	1.302 ± 0.096	1.151 ± 0.091	⋯
2016 Feb 27	1.318 ± 0.122	1.417 ± 0.086	1.376 ± 0.063	1.643 ± 0.133
2016 Mar 26	1.221 ± 0.120	1.441 ± 0.083	1.255 ± 0.064	⋯
2016 Mar 29	1.457 ± 0.228	1.389 ± 0.122	1.468 ± 0.101	1.715 ± 0.141
2016 Apr 03	1.515 ± 0.146	1.447 ± 0.175	1.293 ± 0.089	1.731 ± 0.231
2016 Apr 08	1.570 ± 0.121	1.519 ± 0.133	1.309 ± 0.063	1.445 ± 0.159
2016 May 17	1.255 ± 0.183	1.563 ± 0.100	1.492 ± 0.087	1.998 ± 0.150

Download table as:  ASCIITypeset image

6.1.2. Treatment of Limb Darkening

Our default approach to obtaining the transmission spectrum was to model limb darkening using the quadratic law, with the linear coefficient included as a parameter of the fit and the nonlinear coefficient fixed to a value interpolated from theoretical predictions by Claret (2004). This causes a dependence on theoretical model atmospheres, which is an important consideration for a planet host star as cool as GJ 1132 A. We therefore calculated the planetary radius measurements in the GROND and PISCO bands using an alternate approach: fixing both limb-darkening coefficients.

Figure 5 shows the effect on the results for the GROND bands; the PISCO bands show a similar effect. We find that fixing both limb-darkening coefficients causes the solution to move to higher planetary radius. While the effect is similar for all four bands, it affects r and i more than g and z. It therefore modifies the transmission spectrum at the 1σ–2σ level, in particular the slope seen through the optical wavelength region. As the g and z bands are relatively unaffected, limb darkening cannot explain the anomalously large radius we find in the z band. However, we conclude that observed transmission spectra can be significantly affected by the treatment of limb darkening and urge future studies to check for this possibility in all cases.

6.1.3. Effect of Contaminating Light

Our default solution includes a third-light contribution from a nearby star. This phenomenon can have a significant effect on the transmission spectrum when the spectral energy distribution of the target and contaminant differ (Southworth & Evans 2016), as is the case here. We tested the importance of this effect by calculating solutions for the GROND data while neglecting the third-light contribution.

Once again, the results are represented in Figure 5 and show that the contamination from the nearby star has a negligible effect on the results. The effect is $0.12\sigma$ for the g band and smaller for the other bands. The PISCO data will be similarly (un)affected. We conclude that our transmission spectrum measurements are robust against an amount of contaminating light significantly in excess of that currently known for GJ 1132.

6.1.4. Temporal Variability

We modeled the light curves from all passbands but for each transit individually, in order to check for the presence of temporal variability in the planet. To deal with the complication of having z-band data available for only a subset of the transits, we obtained results for the g, r, and i bands together, and also for all four optical bands when possible. Due to the scatter in the data, we fixed the limb-darkening coefficients for this analysis, using the linear coefficients 0.39 for gri and 0.31 for griz and the quadratic coefficients 0.45 for gri and 0.51 for griz.

Our results are shown in Figure 7. There is no overall trend in measured planet radius during our observations.

Zoom In Zoom Out Reset image size

6.1.5. Stellar Activity

Low-mass stars frequently show dark starspots, which are capable of modifying the transit depth. Starspots occulted by the planet will make the transit shallower, and unocculted starspots will make it deeper (e.g., Ballerini et al. 2012; Oshagh et al. 2013; Tregloan-Reed et al. 2013). GJ 1132 A is a relatively old star and so will not show strong activity, but the fact that a rotation period has been observed (125 days; BT15) means that it does show some starspots. However, there is evidence that many M dwarfs show a large number of small spots (Jackson & Jeffries 2012, 2013), which would greatly reduce the effect of spot activity because the occulted part of the stellar surface would be very similar to the unocculted areas. This is consistent with the lack of starspot features seen in transits of planets in front of M dwarfs (e.g., Mancini et al. 2014; Awiphan et al. 2016).

The very small radius of GJ 1132 b to that of GJ 1132 A means that unocculted starspots would be the dominant contributor to transit depth variations for this system. Because starspots are cooler than the rest of the stellar surface, unocculted spots would have a wavelength-dependent effect on the transit depth. They would cause transits to become gradually deeper as one observed at bluer wavelengths. This does not provide an explanation for the current case, where the z-band radius is significantly larger than the gri-band points. For completeness, we note that unocculted plage could have the opposite effect (see also Oshagh et al. 2014), but this has never been observed. It would also require the plage to have a very clumpy distribution in order to not be occulted during the transits we observed, which is not the situation seen for the Sun, and would not cause an abrupt change in transit depth as found in the z band.

Finally, we note that our observations of GJ 1132 were distributed over 124 days and therefore fortuitously sample the 125-day rotation period of GJ 1132 A very well. Moreover, observations in the different passbands were obtained simultaneously, so temporal variations are unable to affect the relative transit depth measurements. We therefore conclude that our observations are not significantly affected by spot activity in the host star. Having checked for and ruled out issues due to individual transits, contaminating light, and planetary and stellar variability, we conclude that our transmission spectrum is robust.

Our multiband photometric observations allow us to place joint constraints on the interior and atmospheric composition of GJ 1132 b. As discussed above, our observations include measurements of transit depth, and hence the planetary radius, in seven photometric bands. As shown in Figure 4, six of these seven measurements are consistent with a "surface radius" of 1.375 ${R}_{\oplus }$ to within the ∼2σ uncertainties. By "surface radius" we mean the smallest radius that is observationally accessible. However, the measurement in the z band at 0.9 μm differs from a flat spectrum by over 4σ, thereby making it highly statistically significant, and suggests a potential contribution from the planetary atmosphere. Additionally, the z band also overlaps with a strong H₂O absorption band, whereas the remaining six bands probe windows in the H₂O opacity. Therefore, it is possible that these six bands are all measuring the radius of the planetary "surface," while the z band is probing a potentially H₂O-rich atmosphere contributing to a higher measured radius. The importance of using radii measured in different bands to represent the interior versus atmosphere has been suggested previously by Madhusudhan & Redfield (2015), an effect that we are likely witnessing in the present case. Therefore, in what follows, we use a combination of interior and atmospheric models to jointly interpret the data.

The observed mass and radius of GJ 1132 b allow us to investigate the possible interior composition of the planet using internal structure models. The planet mass as shown in Table 4 is 1.63 ± 0.54 ${M}_{\oplus }$. For the radius, we use the baseline value of 1.375 ± 0.16 ${R}_{\oplus }$ discussed above to represent the bulk "gray" radius of the planet. Figure 8 shows model mass–radius curves of homogeneous super-Earths of different compositions spanning Fe, silicates, and H₂O. The models are described in Madhusudhan et al. (2012). We find that the mass and radius are consistent with two broad compositional regimes. First, an exactly Earth-like composition, with 33% iron, 67% silicates, and no volatile layer, is inconsistent with the data within the 1σ uncertainties. However, a composition with a higher silicate-to-iron fraction, including a pure silicate planet, is ostensibly consistent with the data, albeit marginally. On the other hand, the data are also consistent with a large range of H₂O mass fractions between 0% and 100% in our models. In principle, consideration of temperature-dependent internal structure models would lead to larger model radii for the same composition (Thomas & Madhusudhan 2016) and therefore could lower the upper limit on the water mass fraction. Nevertheless, the mass and radius of GJ 1132 b allow for a degenerate set of solutions ranging between a purely silicate bare-rock planet and an ocean planet with a substantial H₂O envelope. The degeneracy between the two scenarios could potentially be resolved using spectroscopic observations of the planetary atmosphere, as discussed below.

Zoom In Zoom Out Reset image size

Considering the transmission spectrum of the planet, we report a tentative inference of H₂O in the planetary atmosphere that provides initial signs of a water-rich world. Figure 9 shows the data and model transmission spectra of GJ 1132 b with different H₂O mixing ratios in an H₂-rich atmosphere; other compositions are explored in the following section. The data are inconsistent with a flat spectrum by over 4σ. As discussed above, the key constraint on the atmospheric composition of the planet is governed by the z-band measurement, which shows a substantially higher transit depth relative to the baseline. We model the atmospheric transmission spectrum of the planet using the exoplanetary atmospheric modeling method of Madhusudhan & Seager (2009) and Madhusudhan (2012). Given the limited number of data points, we did not embark on a full retrieval exercise, but instead systematically investigated a grid of model compositions. In this section, we consider models comprising only gaseous H₂ and H₂O at different mixing ratios to explore the potential contribution of H₂O to the z-band measurement. We nominally fixed the temperature structure to be isothermal at 600 K, representative of the equilibrium temperature of the planet; nominal variations in the temperature do not change our conclusions significantly. We find that the best model fits to the data, particularly the z-band point, are obtained for H₂O volume mixing ratios of 1%–10%, as shown in Figure 9. On the other hand, the remaining data that are consistent with a flat spectrum are fit relatively easily for a wide range of models, as the corresponding photometric bands largely probe windows in opacity.

Zoom In Zoom Out Reset image size

The amplitude of the observed z-band feature is physically plausible for a range of compositions. The difference between the z-band radius and the continuum provides an estimate of the thickness of the observable atmosphere (${H}_{\mathrm{atm}}$) of 0.22 ± 0.06 ${R}_{\oplus }$ or $\sim 1400\pm 400$ km. Considering an isothermal atmosphere at the equilibrium temperature (644 K), the atmospheric scale height (${H}_{\mathrm{sc}}$) for an H₂-dominated atmosphere (mean molecular mass $\mu =2$) is ∼300 km. The value of ${H}_{\mathrm{sc}}$ for 10% H₂O ($\mu =3.6$) and 20% H₂O ($\mu =5.2$) is 167 and 115 km, respectively, assuming that H₂ occupies the remaining fraction. Therefore, the atmospheric height can be explained by ∼5 scale heights of an H₂-dominated atmosphere given a strong absorber in the z band. Similarly, ${H}_{\mathrm{atm}}$ can also be explained by $\lesssim 10$ scale heights of a 10% H₂O atmosphere. Conversely, considering ∼8–10 scale heights expected for a saturated spectral feature (Madhusudhan & Redfield 2015) at the same temperature, the mean molecular mass of the atmosphere is constrained to be $\mu \sim 2.8$–5.5. Such a μ is possible in an H₂-rich atmosphere with ∼10%–20% H₂O or corresponding fractions of other heavy molecules, again assuming that a strong absorber in the z band is present in the atmosphere. Future observations with the Hubble Space Telescope (HST) and James Webb Space Telescope in the near-infrared could further constrain the presence and composition of such an absorber.

The combined interior and atmosphere models fitting the data are consistent with a water-rich envelope in GJ 1132 b. And given that the data are inconsistent with a flat spectrum, the atmosphere is unlikely to be of very high mean molecular mass, e.g., 100% H₂O or CO₂, or one with thick high-altitude clouds. However, since the inference is based primarily on one photometric data point (albeit obtained from light curves of six transits), future observations are critical to validate our current finding. We predict that a transmission spectrum of GJ 1132 b observed with HST/WFC3 in the 1.1–1.8 μm spectral range should be able to detect the water absorption in its atmosphere. On the contrary, if such an HST spectrum does not detect an H₂O feature, then it might suggest the presence of another molecule in the atmosphere that we have not accounted for in our present models, or the possibility of time-variable events in the planetary atmosphere. It is also necessary that additional observations of the planetary transit in the z band be conducted to bolster the current finding.

7.1. Transmission Spectra Calculations Using PetitCODE

Using the newest version of the petitCODE (Mollière et al. 2015, 2016), we performed a second, independent exploration of the atmospheric properties of the planet. petitCODE self-consistently calculates the radiative-convective equilibrium structures of irradiated or self-luminous exoplanet atmospheres. Molecular and atomic line opacities, cloud opacities, and H₂–H₂ and H₂–He collision-induced absorption (CIA) are taken into account. The code also includes molecular Rayleigh and cloud particle scattering.

For these calculations we used the planetary parameters as defined in Table 4. We assumed a stellar effective temperature of 3270 K from BT15. We calculated our standard suite of models as defined in Mollière et al. (2016) for this planet, consisting of one fiducial model without clouds, two models at half and twice the solar C/O ratio, respectively, and nine different cloud model approaches as defined in Table 2 of Mollière et al. (2016). For all models (three clear and nine cloudy) we calculated atmospheric structures and spectra at four different scaled-solar chemical compositions, corresponding to atmospheric metallicities of $\left[\tfrac{\mathrm{Fe}}{{\rm{H}}}\right]=0$, 1, 2, and 3.

We considered Na₂S and KCl clouds only, because higher-temperature cloud species can most likely not be mixed up to the higher layers of the atmosphere (Charnay et al. 2015; Parmentier et al. 2016). Cloud models 1–4 represent models of varying cloud thickness using the model by Ackerman & Marley (2001) and adopting ${f}_{\mathrm{sed}}$ values of 0.01–3, where ${f}_{\mathrm{sed}}$ is the ratio of the particle-mass-averaged settling speed and the atmospheric mixing velocity. Cloud models 5–9 represent extended clouds with monodispersed size distributions. The particle size is fixed at 0.08 μm, which leads to Rayleigh scattering of the clouds in the optical. The cloud mass fraction is equal to the mass fraction derived from equilibrium chemistry, but not larger than a maximum value ${X}_{\max },$ which decreases from model 5 to 7. ${X}_{\max }$ can thus be thought of as a simple parameterization of the settling strength. Cloud models 8 and 9 adopt the same parameter choice as cloud model 6, but additionally include iron clouds (model 8) or a spherical, homogeneous cloud particle shape (model 9). The cloud models therefore cover the parameter space from larger-particle clouds (models 1–4), leading to flatter transmission spectra in the optical, to small-particle clouds (models 5–9), which lead to Rayleigh scattering. Model 8 will not exhibit Rayleigh scattering if iron can condense within the atmosphere, due to the strong optical absorption of iron.

Each of the model spectra was compared to the observed transmission spectrum and the level of agreement determined. This was done by identifying the lowest ${\chi }^{2}$ value between the observed transmission spectra and passband-integrated model values from petitCODE, while allowing for an overall shift in radius because the spectra were calculated assuming a planetary base radius of 1.43 ${R}_{\oplus }$ at 10 bar pressure. Faced with a multitude of choices over which data to consider, we defaulted to using the GROND griz points with one limb darkening coefficient fitted and the GROND JHK points with fixed limb darkening coefficients. This is because there is a disagreement between the planet radius in the GROND and PISCO g bands, and the former was obtained from nine transits versus the single transit for the latter. We prefer results obtained with fitted limb darkening coefficients because of limitations in the theoretical understanding of the atmospheres of stars as cool as GJ 1132 A. Alternative choices of data points will be discussed below. Selected models are shown in Figure 10.

Zoom In Zoom Out Reset image size

We found that the best fit to the GROND grizJHK transmission spectrum is obtained for the petitCODE model with $\left[\tfrac{\mathrm{Fe}}{{\rm{H}}}\right]=2$ and a twice-solar or solar C/O ratio (${\chi }^{2}=18.6$ and 19.1, respectively, for seven data points), although acceptable agreement also occurs for the models with $\left[\tfrac{\mathrm{Fe}}{{\rm{H}}}\right]=1$ and $\left[\tfrac{\mathrm{Fe}}{{\rm{H}}}\right]=3$ (${\chi }^{2}$ from 19.9 to 20.7 depending on C/O ratio). The models with $\left[\tfrac{\mathrm{Fe}}{{\rm{H}}}\right]=3$ have muted spectral features, because the high mean molecular mass leads to a low atmospheric scale height. They match the H-band radius better, at the expense of a poorer agreement with the z-band radius. A featureless spectrum (i.e., a straight line) is a worse fit with ${\chi }^{2}=22.1$. The ${\chi }^{2}$ values are all quite large and reflect the difficulty of finding a good agreement between the observed and theoretical spectra. Including clouds in the models serves to increase the ${\chi }^{2}$ because the flattened spectral features match the observed transmission spectrum less well.

If we restrict the analysis to the GROND griz data points, as these are the most reliable measurements, the models with $\left[\tfrac{\mathrm{Fe}}{{\rm{H}}}\right]=1$ and $\left[\tfrac{\mathrm{Fe}}{{\rm{H}}}\right]=2$ are the best matches. The greatest agreements are for $\left[\tfrac{\mathrm{Fe}}{{\rm{H}}}\right]=2$ and a twice-solar C/O ratio (${\chi }^{2}=8.9$ for four data points), followed by $\left[\tfrac{\mathrm{Fe}}{{\rm{H}}}\right]=1$ and a solar (${\chi }^{2}=10.0$) or twice-solar (${\chi }^{2}=10.5$) C/O ratio. The straight-line fit is clearly disfavored with ${\chi }^{2}=16.9$. The slight preference for the twice-solar C/O ratio is because a higher C/O ratio enhances the amount of CH₄, an effect that outweighs the loss of H₂O and hence leads to a larger radius in the z band. Even these models underpredict the detected variation of planet radius with wavelength, which also means that a cloudy or featureless spectrum is strongly disfavored. We therefore conclude that our results are best explained by a clear atmosphere with $\left[\tfrac{\mathrm{Fe}}{{\rm{H}}}\right]=2$ and strong z-band opacity contributed by CH₄ or H₂O.

Turning now to the transmission spectrum from the GROND data with all limb darkening coefficients fixed, we find that the best models have $\left[\tfrac{\mathrm{Fe}}{{\rm{H}}}\right]=2$, with ${\chi }^{2}=19.4$ (seven data points) for a solar C/O ratio, 19.5 for twice-solar C/O, and 19.8 for half-solar C/O. $\left[\tfrac{\mathrm{Fe}}{{\rm{H}}}\right]=3$ with twice-solar C/O is nearly as good (${\chi }^{2}=20.2$), as is $\left[\tfrac{\mathrm{Fe}}{{\rm{H}}}\right]=1$ with half-solar C/O (${\chi }^{2}=20.8$). A straight-line fit has ${\chi }^{2}=22.1$. The quality of our data and limitations in the limb darkening treatment means that we are able to infer that the atmosphere of GJ 1132 b is best represented by models with $\left[\tfrac{\mathrm{Fe}}{{\rm{H}}}\right]=2$, but the C/O ratio is not usefully constrained.

GJ 1132 is a benchmark nearby system containing a low-mass planet transiting a late M dwarf. We have presented extensive photometry of the system comprising light curves of nine transits observed simultaneously in the griz optical and JHK near-IR passbands. We have analyzed these and literature data to determine the physical properties of the system. We find that the planet is larger than previously thought, 1.43 ± 0.16 ${R}_{\oplus }$ versus 1.16 ± 0.11 ${R}_{\oplus }$, from a methodological approach that relies more on observations of the system and less on empirical calibrations of the properties of low-mass stars. The planet's measured mass and radius are consistent within 1σ with theoretical predictions for a planet composed of silicates or water; a 100% iron composition gives a radius too small by ∼2σ.

Our repeat observations allowed us to check for variability in the measured planet radius, and two of the transits do indeed yield radii that are modestly discrepant with measurements from other transits. This could indicate excess scatter among the results, starspots on the stellar surface, unidentified systematic effects in our data, or the presence of variability in the planet's atmosphere (e.g., Armstrong et al. 2016).

We have constructed an optical–infrared transmission spectrum of GJ 1132 b by modeling all light curves with a consistent geometry. We find an increased planet radius in the z band, to a significance level of 4σ, indicative of atmospheric opacity due to water, methane, or another unidentified source. Detailed investigation of the resulting error bars was enabled by the observation of nine transits. We find that our results are robust against the rejection of individual transits or the inclusion of contaminating light from a nearby star. The treatment of limb darkening is more concerning, as it affects the results in the r and i bands at the level of 2.7σ and 3.5σ, respectively. We urge fellow researchers to consider this issue in similar analyses, especially for very cool stars where theoretical limb darkening coefficients are less reliable.

The transmission spectrum was modeled using the atmospheric models of Madhusudhan & Seager (2009), with the finding that H₂O likely causes the enlarged z-band radius of the planet. The best fits to the observations are found for H₂O volume mixing ratios of 1%–10%, implying a water-rich atmospheric composition that would cause observable spectral features in a 1.1–1.8 μm transmission spectrum obtained using HST/WFC3. From simulations of the atmosphere of GJ 1132 b, Schaefer et al. (2016) found that the presence of H₂O implied either an H₂ envelope or low UV flux from the host star early in the lifetime of the system, and the ongoing presence of a magma ocean on the planet's surface.

We also calculated theoretical spectra using petitCODE (Mollière et al. 2015), which yield similar results except for the finding that the large z-band radius is explicable by an enhanced abundance of CH₄. A high metallicity of $\left[\tfrac{\mathrm{Fe}}{{\rm{H}}}\right]=2$ is preferred, depending on the data points considered, which is in line with the mass–metallicity correlation seen for more massive planets (Kreidberg et al. 2014a; Mordasini et al. 2016). A straight line is a much poorer fit to the transmission spectrum, confirming that we have detected the atmosphere of a 1.6 ${M}_{\oplus }$ planet.

We advocate extensive further observations to refine and extend our understanding of the GJ 1132 system. High-precision optical light curves from large telescopes would be able to confirm or disprove the larger radius of the planet in the z band and shed light on the discrepancy seen in the g band. Intermediate-band photometry at 900 nm or bluer than 500 nm would enable finer distinctions to be made between competing model spectra and a clearer understanding of the chemical composition of the planetary atmosphere. The planet's mean density measurement is also hindered by the weak detection of the velocity motion of the host star, an issue that could be ameliorated with further radial velocity measurements using large telescopes. Finally, infrared transit photometry and spectroscopy should allow the detection of a range of molecules via the absorption features they imprint on the spectrum of the planet's atmosphere as backlit by its host star.

Our results show that a 1.6 ${M}_{\oplus }$ planet with an equilibrium temperature of 650 K is capable of retaining an extensive atmosphere. The atmosphere contains multiple molecular species and has likely persisted for many gigayears since the formation of the system.

J.S. acknowledges financial support from the Leverhulme Trust in the form of a Philip Leverhulme Prize. We thank BT15 for making their photometric data available. The following Internet-based resources were used in research for this paper: the ESO Digitized Sky Survey; the NASA Astrophysics Data System; the SIMBAD database and VizieR catalog access tool operated at CDS, Strasbourg, France; and the arXiv scientific paper preprint service operated by Cornell University.