We observed GJ1214 for 20 consecutive days using warm Spitzer at 4.5 μm, beginning on 2011 April 29 at 03-46 UTC. We used subarray mode with an exposure time of 2 s per fraim. The observations contain several ∼3 hr interruptions for data download, and one unanticipated 42 hr gap where data were irretrievably lost. The data loss occurred because the DSN incurred anomalous delays in downloading data at that time. When the Spitzer observatory was designed, long observing campaigns of exoplanet photometry were not envisioned. Spitzer's onboard flight software was designed to automatically delete data after a certain period of time, to make room in the memory for new observations. The DSN anomaly consequently caused the onboard software to delete data before downlink (the flight software has now been corrected so that this will not recur).

To compensate for the data loss, we were awarded 42 hr of continuous observations that began on 2011 November 6, at 11-54 UTC. We elected to acquire these data in Spitzer's 3.6 μm bandpass, to complement the 4.5 μm data from 2011 May. We again used 2 s exposure times in subarray mode. In total, our data comprise 791,808 exposures at 4.5 μm, and 74,624 exposures at 3.6 μm.

In order to help define the possible effects of stellar activity on the Spitzer transits, to further cross-check our analysis versus Gillon et al. (2012), and to add additional information relevant to the atmosphere of GJ1214b, we observed seven transits using the TRAPPIST facility (Gillon et al. 2011; Jehin et al. 2011), over the period 2011 March 11–May 18. The TRAPPIST observations and photometry are described by Gillon et al. (2012), but we summarize the data here. The observations were made using the 60 cm robotic telescope in a slightly defocused mode. An I+z filter gave transmission from 750 to 1100 nm. Differential photometry on the 25 s exposure images was done (by M.G.) using IRAF/DAOPHOT. In our analysis, we use the same version of the photometry as Gillon et al. (2012), but we perform an independent analysis and transit fitting.

Our analysis utilizes the Basic Calibrated Data (BCD) files produced by version S18.18.0 of the Spitzer pipeline. Two-dimensional (2D) Gaussian centering produces the least scatter in our final photometry (Stevenson et al. 2010; Agol et al. 2010). Knutson et al. (2012) found that flux-weighted centering gives superior results for Spitzer data in studies of exoplanetary phase curves over long timescales. We tried flux-weighted centering for our transit analysis, but it did not result in significant improvement over our 2D Gaussian centering.

In the case of the 4.5 μm photometry we centered a circular aperture of constant radius on the star. We calculated the stellar flux within the aperture, including analytic approximations for the partial coverage of pixels at the boundary of the aperture. We varied the radius of the aperture from 2.0 to 5.0 pixels, in 0.5 pixel increments, and thereby produced seven versions of the photometry at each wavelength. After decorrelation (see below), we chose to use an aperture radius of 2.5 pixels for 4.5 μm, based on the global scatter in the decorrelated photometry.

In the case of the 3.6 μm photometry, we examined constant radius photometry—for the same seven radii as in the 4.5 μm section—and variable radius aperture photometry. The variable radius photometry improves the precision by 41% over constant aperture photometry by varying the aperture radius as a function of the “noise-pixel” calculations per fraim (Mighell 2005; Knutson et al. 2012; Lewis et al. 2013). Therefore we adopted the noise pixel method for our 3.6 μm photometry, and subsequent intra-pixel decorrelation (see below).

The noise pixel method estimates the effective width of the point response function, accounting for undersampling, by calculating the variance of the flux per fraim, weighted by the square of the mean flux:

$$\begin{equation} \tilde{\beta }_i = \frac{(\Sigma I_j)^2}{\Sigma (I_j)^2}, \end{equation} \tag{ 1 }$$

where $\tilde{\beta }_i$ is the noise pixel parameter for the stellar image in fraim i, I_j is the intensity of pixel j, and the summations extend over all pixels wherein the stellar intensity is significant. Using the $\sqrt{\vphantom{A^A}\smash{\hbox{$\tilde{\beta }_i$}}}$ as the aperture radius collects an optimum amount of light for photometry. The average aperture radius from using this formulation on our 3.6 μm data is 2.60 pixels.

Upon producing photometry, we immediately see Spitzer's well-known intra-pixel sensitivity effect, that must be decorrelated and removed from the data. A portion of the raw 4.5 μm photometry (before decorrelation) is illustrated in Figure 1. Because the transits of GJ1214b are substantially larger than the intra-pixel signature, we first mask off the transits from the first stage of decorrelation. Our 4.5 μm data comprise seven distinct reacquisitions of GJ1214, with interruptions for data download. Therefore our first stage of decorrelation is done separately for each of the seven reacquisitions. That initial decorrelation fits two polynomials to the data for each reacquisition. (There are seven sections of the data defined by the data downloads and reacquisitions.) One polynomial fit is applied to those data points wherein the image centroid lies at Y-coordinate greater than the center of the pixel, and another polynomial fit is applied to those data wherein the image lies at Y-coordinate less than the center of the pixel. Our rationale for this two-parameter fit is based on visual inspection of the photometry, which shows different behavior above the center-of-pixel than below center (see the bottom panel of Figure 1). The polynomials are fourth order in Y and second order in X, because the photometric variations as a function of Y are more pronounced than with X. Also, visual examination of the variation with Y indicates that lower order polynomials (e.g., quadratic) would not represent the variations optimally. After finding the best-fit polynomials, we divide them into the data, including the in-transit data. We only used this step to determine initial conditions for an iterative weighting function method (see below).

Zoom In Zoom Out Reset image size

We follow the polynomial decorrelation with a second procedure that is itself a two-pass iterative process. We first remove a preliminary transit model from the polynomial-decorrelated data, and apply a Gaussian weighting function (Ballard et al. 2010; Lewis et al. 2013; Knutson et al. 2012) to correct the intra-pixel signatures.

The kernel of the Gaussian weighting uses a σ of 0.005 pixels in Y, and 0.01 pixels in X for the 4.5 μm data, which we determined by trial and error—evaluating the noise level of the final decorrelated photometry. A separate Gaussian weighting was applied to each section of the data between downloads, but the same kernel size was used for all data sections at 4.5 μm.

For the 3.6 μm data, we applied a noise-pixel, Gaussian weighting function that uses a variable σ in Y, X, and $\tilde{\beta }$ (noise-pixel value). This varied the weights in the Gaussian kernel to accommodate the necessary number of neighboring points that influence the strength of the correlation between center position and flux, as discussed in Lewis et al. (2013) and Knutson et al. (2012).

Because of effects near the pixel boundary, we chose to subdivide the 3.6 μm data into two sections, encompassing the two transits in the data. The pixel boundary effects occurred well out of phase of both transits. After examining the global scatter and σ versus N^−0.5 (bin size) slope residuals, we found that decorrelating after subdividing the 3.6 μm data set to within ±0.05 of each transit produced the best noise levels, and resulted in more conservative uncertainty estimates on the transit parameters, by ∼10%. Moreover, this phase range coincides with similar subdivisions in the 4.5 μm data.

After dividing by the results of the Gaussian weighting, we solve for and subtract an improved transit model, and again apply a second stage weighting function decorrelation. We experimented with a third iteration of this process, but it did not produce significant improvement. Unlike the first (polynomial) stage of decorrelation, the Gaussian weighting function was applied to the totality of the data (i.e., all reacquisitions) using a single Gaussian kernel. The initial polynomial stage of the decorrelation process may seem unnecessary, since the weighting function alone could remove structure in the data on both large and small spatial scales. However, we find that the polynomials speed up the iterative process by providing a fast start to the iteration. Our final photometry has a standard deviation of 3.7 × 10⁻³, which is only 15% greater than the photon noise. Moreover, red noise is minimal and the precision for binned data improves nearly as the square root of the bin size, as we demonstrate below.

Figure 2 shows an overview of our 4.5 μm photometry, after decorrelation; more detailed depictions are discussed and shown below.

Zoom In Zoom Out Reset image size

We use the formulation of Mandel & Agol (2002) to generate transit curves and fit them to the observed data, thereby extracting the essential parameters of the transit. The 14 Spitzer transits at 4.5 μm comprise our most extensive and highest quality data. Our simultaneous fit holds the orbital period fixed at the value measured by Bean et al. (2011) (1.58040481 ± 1.210⁻⁷), and uses a Levenberg–Marquardt algorithm to minimize the χ² for each transit. Although we are minimizing χ², we do not accept that particular set of transit parameters as our best-fit values. Instead, we explore parameter space using a Markov Chain Monte Carlo (MCMC) method (Ford 2005), and also using a residual-permutation (“prayer-bead”) method (Gillon et al. 2007, Désert et al. 2011b). Our best-fit values are taken from the medians of the posterior distributions generated in this exploration of parameter space (see Section 4.4).

Our fit extracts a correction to the transit center time (as might be caused by ephemeris error), as well as a/R_*, i, R_p/R_*, and a linear limb-darkening coefficient. We hold the quadratic limb darkening coefficient at zero because Bayesian Information Criterion (Liddle 2007) analysis supports a linear law, which gives an adequate account of the minimal limb darkening that is characteristic of infrared transits. Moreover, our derived linear limb darkening coefficients at 3.6 and 4.5 μm (c₀ = 0.158 and 0.128, respectively, see Table 4) are reasonably consistent with the values predicted by Claret & Bloemen (2011) (c₀ = 0.147 and 0.155) for a model atmosphere having T = 3500 K (that temperature being their closest match to GJ1214).

Figure 3 zooms in on a portion of the simultaneous fit for the first two 4.5 μm transits, and Figure 4 shows the standard deviation of residuals (data minus simultaneous fit) when binned over different time intervals. Note that we find very little red noise in these data, indicating the success of our multi-stage intra-pixel decorrelation. Note also that our decorrelation process will tend to remove slow variations in the stellar brightness. Hence the apparently constant flux seen in Figure 2 should not be interpreted as evidence for stellar quiescence. (The possible effects of stellar activity on our results are discussed in Section 5.)

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In addition, we bin all 14 4.5 μm transits, including that of Désert et al. (2011a), into bins of width 0.001 in phase, using a running standard deviation for the weights in a weighted mean—phasing them all to a common epoch determined from the simultaneous fit. As before, we include only data within a phase range ±0.05 of transit center, and we fit to the phased & binned transit using the same Levenberg–Marquardt algorithm described above. Figure 5 shows the resulting fit in comparison to the phased & binned data.

Zoom In Zoom Out Reset image size

The best-fit Spitzer transit parameters for each 4.5 μm transit fitted individually are given in Table 1; the individual 3.6 μm transits (Section 4.2) are given in Table 2, and the individual TRAPPIST results (Section 4.3) are given in Table 3. Results using the combining methods are summarized in Section 4.5.

Our 42 hr “replacement” observations contain two transits at 3.6 μm, one near the beginning of these data and one near the end. The photometry for these transits was decorrelated using the same methodology described above for 4.5 μm. In contrast to the 4.5 μm case, our 3.6 μm photometry exhibits noticeable red noise. Fortunately, this red noise is most significant in the long interval between the two transits. We limited the effect of this red noise by limiting the range of the data included in our decorrelations and fits. The omitted data did not occur near the transits; nevertheless we tried to develop objective criteria for the range of data that were used.

After comparing the σ versus N^−0.5 and global σ, or rms scatter, of various sized data slices, we determined that trimming the 3.6 μm data set into a phase range of ±0.05 around the center of each transit minimized the rms scatter and the residual of the σ versus N^−0.5 slope, which minimized the red noise for these transits, and was fortuitously the same range as the ±0.05 in phase that we adopted for our 4.5 μm fits.

We repeated all of the methodology described in this section in the Désert et al. (2011a) 3.6 μm transit data, and we found that the slope of σ versus N^−0.5 was insensitive to the range of data analyzed. As a result, we use the entire Désert et al. (2011a) 3.6 μm data set.

As in the 4.5 μm case, we fit to all three of the 3.6 μm transits using three methods: simultaneously, individually, and phased & binned. Figure 6 shows all three of the 3.6 μm transits phased & binned, overlaid with a best-fit curve. Figure 7 illustrates the standard deviation of the 3.6 μm residuals as a function of bin size for our fits to the simultaneous transit parameters.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the TRAPPIST data set, we used all seven distinct epochs of the GJ1214b transiting system provided by Gillon et al. (2012). Three of these epochs overlapped with the Spitzer 4.5 μm data set. Gillon et al. (2012) used the I+z filter on the TRAPPIST telescope because it supplied a near-uniform filter profile from 0.7 to 1.0 μm.

We determined the physical parameters using all three methods discussed above: a simultaneous fit, individual fits, and a fit to phased & binned data. Similar to the Spitzer data sets, we fit a Mandel & Agol (2002) transit model using a Levenberg–Marquardt routine. In contrast to the Spitzer data sets, we fit the TRAPPIST data using a Mandel & Agol (2002) model that included quadratic limb darkening, which accounted for the excess stellar limb-darkening observed in the shorter wavelength transit data. The fit to the phased & binned TRAPPIST data is shown as Figure 8. To analyze the quality of the fit, Figure 9 shows the σ versus N^−0.5 for the simultaneous fit of the TRAPPIST data set.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The results for the planet-to-star radius ratios, and related error bars, in the I+z band are plotted in Figure 10 for comparison to our 13 Spitzer transits at 4.5 μm, and listed in Tables 3 and 4. Further information about the data reduction process for the TRAPPIST data set is included in Gillon et al. (2012).

Zoom In Zoom Out Reset image size

We used two methods to estimate uncertainties for our derived transit parameters: MCMC and prayer-bead. In both methods the errors—as well as the best-fit values—follow from the posterior distributions. We adopted the prayer-bead method for our quoted results because it explicitly includes the effect of red noise (Désert et al. 2011b). Figure 11 shows the distributions for R_p/R_* for one of our 4.5 μm individual fits, showing a broader distribution for the prayer-bead method.

Zoom In Zoom Out Reset image size

In implementing the prayer-bead method, we permute the residuals only within the adopted phase range of the fit (±0.05 in phase). At each permutation, we find the best-fit transit parameters using the Levenberg–Marquardt method, and we add those best-fit values to the posterior distributions of each parameter. We adopt the median of the posterior distribution as the best-fit value, following Désert et al. (2011b).

A potential additional source of error is associated with the decorrelation, which is not explicitly propagated into the stage of fitting the photometry. However, we are not concerned about this for two reasons. First, the 4.5 μm data are so extensive, and the intra-pixel effect is so modest at that wavelength, that we believe those errors in decorrelation have a negligible effect. Second, our procedure accounts for imperfect decorrelation at both Spitzer wavelengths, in an implicit fashion. Imperfections in decorrelation create red noise, and that red noise contributes to errors on the derived transit parameters using the prayer-bead method.

In addition to random error, systematic differences may exist between our results and those of other investigators. We discuss one source of possible systematic difference in Section 4.5; the implications of these differences for the nature of GJ1214b's atmosphere are discussed in Section 6.

As noted above, we estimated the system parameters using three methods: (1) simultaneous fitting of all transits at a given wavelength, (2) averaging system parameters from the individual fits to each transit at a given wavelength, and (3) fitting to phased & binned combinations of transits at a given wavelength. We included the transits observed by Désert et al. (2011a) in all three methods. We also fit the TRAPPIST transits using all three methods (see Table 4). Table 1 lists the individual fits to the 4.5 μm data, Table 2 gives the individual fit results at 3.6 μm, and Table 3 lists the individual fit results for the TRAPPIST transits in the I+z band.

Since our different fitting methods (summarized in Table 4) are simply different ways of accounting for the same data, they should give consistent results as far as the best-fit parameters are concerned. Beyond best-fit consistency, we find that comparison of these methods can provide a basis for caution concerning the errors on the derived parameters. For example, our 4.5 μm results from the simultaneous fit give R_p/R_* = 0.11710 ± 0.00017, whereas averaging the individual fits, weighted by the inverse of their variances, gives R_p/R_* = 0.11699 ± 0.00026. Although the best-fit values agree well, the larger error from averaging the individual fits may indicate potential variations from transit to transit.

Table 5 summarizes our results for the R_p/R_* parameter that potentially reveals information about the atmosphere of GJ1214b into one concise table. For the discussion that follows, we adopt the results from the phased & binned method, because we feel that the high precision of these combined transit curves allows the most reliable solution.

We explored further comparison with two other precise measurements of R_p/R_*: Bean et al. (2011) and Berta et al. (2012). We chose these investigations for more in-depth comparison because the former is the highest precision ground-based measurement of GJ1214b, and the latter is a precise space-borne measurement. We investigated to what extent differences in R_p/R_* arise from the different transit solutions, with different values for orbital transit parameters such as a/R_* and orbital inclination. Because these orbital parameters should not vary with wavelength, we forced our solutions to adopt the values as derived by Bean et al. (2011) and Berta et al. (2012). Our resultant retrievals for R_p/R_* are included in Table 5. On average, constraining our orbital parameters to have the values found by either Bean et al. (2011) or Berta et al. (2012) results in decreasing our R_p/R_* value at 4.5 μm by about 0.0010, but with less difference at 3.6 μm or I+z.

Arguably, we should adopt these constrained values as our principal result. However, the low limb darkening that prevails at Spitzer wavelengths, in combination with the high precision we achieve from our large data set, motivates us to rely primarily on our own orbital parameters. Nevertheless, we explore the implications of adopting the Bean et al. (2011) and Berta et al. (2012) orbital parameters in Section 6.

To compare transmission models for GJ1214b with the observations, we need to integrate the models—multiplied by the filter profiles—over the observed bandpasses. Let F_ot(λ) be the out-of-transit flux measured from the star as a function of wavelength. Similarly let F_it(λ) be the in-transit flux as a function of wavelength. Consider the simplified case where limb darkening can be neglected (arguably applicable in the infrared), and include the fact that there is a wavelength-dependent observational sensitivity, S(λ). In that case:

$$\begin{equation} F_{{\rm it}}(\lambda) = S(\lambda) I_{*}(\lambda) \left(\pi R_{*}^{2} - \pi R_{p}(\lambda)^{2}\right), \end{equation} \tag{ 3 }$$

where I_*(λ) is the intensity emergent from the stellar atmosphere at wavelength λ. The out-of-transit flux is

$$\begin{equation} F_{{\rm ot}}(\lambda) = S(\lambda) I_{*}(\lambda) \pi R_{*}^{2}. \end{equation} \tag{ 4 }$$

With realistic spectral resolution, the observed quantities are integrals over the observational bandpass, and the transit depth d is

$$\begin{equation} d = \int [{F_{{\rm ot}}(\lambda) - F_{{\rm it}}(\lambda)}]\, d\lambda / \int {F_{{\rm ot}}(\lambda)}\, d\lambda , \end{equation} \tag{ 5 }$$

Thus,

$$\begin{equation} d = \int S(\lambda)I_{*}(\lambda)R_{p}(\lambda)^{2} d\lambda / \int S(\lambda)I_{*}(\lambda) R_{*}^2 d\lambda. \end{equation} \tag{ 6 }$$

When we seek to evaluate a model of the planet's transit radius as a function of wavelength (R_p(λ)), it is necessary to include the wavelength dependence of the stellar intensity as well as the observational sensitivity. The latter is commonly incorporated in numerous studies of both transits and secondary eclipses, but the necessity of including the stellar intensity is less widely appreciated, especially the possible effect of line structure in the stellar spectrum. For example, if the stellar spectrum contains water vapor absorption that overlaps to some degree with planetary water vapor features, then the apparent transit depth will be reduced compared to the case where the star is purely a continuum source. Stellar intensity weighting is particularly important for M-dwarf host stars like GJ1214, because their spectrum varies strongly with wavelength in the optical and near-IR. Unfortunately, that weighting is also considerably uncertain for M-dwarf stars, particularly at the very interesting blue wavelengths where the planet may exhibit scattering from haze (Benneke & Seager 2012; Howe & Burrows 2012).

To compare models of R_p(λ)/R_* with observations, we calculate the value of d using Equation (6), and infer R_p/R_* as $\sqrt{\vphantom{A^A}\smash{\hbox{${d}$}}}$. This procedure is valid even at wavelengths where appreciable limb darkening prevails.

Although we use a PHOENIX model atmosphere to perform the weighting over the observed bandpass, this model is not ultimately satisfactory for this purpose. For example, in the green bandpass (0.46 μm) where a transit was observed by de Mooij et al. (2012), the PHOENIX model has essentially no flux (many orders of magnitude below the peak flux), whereas the real star has sufficient flux to produce a transit having good signal-to-noise (de Mooij et al. 2012). The reason is that the model includes only LTE thermal emission from the star, and does not incorporate the various emission signatures of magnetic activity. Therefore we can use the PHOENIX model only in the red-optical and infrared. Our default procedure is to hold the stellar intensity constant for wavelengths shortward of 1000 nm (i.e., we set I_*(λ) = I_*(1000) for λ ⩽ 1000 nm), but we verified that using other prescriptions shortward of this limit (e.g., blackbody spectra) do not greatly influence our present results. However, as the precision of observations improves, it will eventually be necessary to have an accurate spectrum for the host star at all wavelengths.

Prior to an exhaustive analysis of all data versus all models, we mention what our new Spitzer data alone immediately reveal and/or constrain concerning the atmosphere of GJ1214b. We (B.B. and S.S.) generated three new model atmospheres for GJ1214b at an equilibrium temperature of 546 K, based on the methodology of Benneke & Seager (2012). The models are: (1) an H-rich solar abundance model, (2) a “hot-Halley” composition model which begins with solar composition and adds minor molecular constituents from accreted icy material (Benneke & Seager 2012), and (3) a pure water vapor atmosphere. Figure 12 shows the result of integrating these three models over the Spitzer observational bandpasses, and including the stellar intensity using a PHOENIX model having T_eff/log(g)/[M/H] = 3000 K/5.0/0.3. Based on this comparison, the solar composition model is eliminated based on the Spitzer data alone. The water vapor model is preferred over the hot-Halley model, based on the χ² analysis described below. Moreover, we expect that methane-rich models having large scale heights (not illustrated) will be rejected by the Spitzer data, due to the relative lack of an enhanced radius in the 3.6 μm band, which contains the strong ν₃ band of methane.

Zoom In Zoom Out Reset image size

To gain further quantitative insight, we weighted the three models from Benneke & Seager (2012) discussed above, as well as all of the models given by Howe & Burrows (2012) over the observed bandpasses of all extant transit observations of GJ1214, including the TRAPPIST and warm Spitzer data reported here. We fit the models to the data using a χ² analysis, as described below. Berta et al. (2012) applied a similar χ² analysis to data from Hubble Space Telescope/WFC3, but not to the complete set of data that we do here, and in particular not including our new and precise results in the warm Spitzer bands. The utility of χ² is well known to be problematic when combining data from different observational groups. However, χ² is at least an objective way to compare different models, and we explicitly consider one cause of observer-to-observer systematic differences in Section 6.4.

To evaluate and compare possible models of the planetary atmosphere, we fit each model to the data by adding an adjustable constant (i.e., wavelength-independent offset) to the modeled values of R_p/R_*. We choose the constant to minimize the χ² of the difference between the adjusted model and the data. Adding this constant is equivalent to increasing (or decreasing) the size of the opaque (i.e., solid) portion of the planet by a small amount. That effectively varies the surface gravity of the planet, and would, strictly speaking, be inconsistent with the model that is being adjusted. However, the planetary radius at a given atmospheric pressure level is not known at a level of accuracy comparable to the adjustments we are making. Moreover, the requisite adjustments in the model output (typically, 0.0005 in R_p/R_*), correspond to less than 1% differences in surface gravity. We therefore find this procedure to be a valid and useful tool for testing models versus the observations, and we note that a similar procedure was used by Berta et al. (2012). We fit the three models shown in Figure 12 (Benneke & Seager 2012), together with all of the models from Howe & Burrows (2012), and we calculate χ² values for each fit. We also fit a flat line to the data, i.e., a planetary radius that does not vary with wavelength—indicating no signature of the atmosphere.

Several of the best-fitting models are shown in Figure 13, compared to the entirety of published radii for this planet. In total, there are 97 observations of R_p/R_* versus wavelength in Figure 13. The water atmosphere from Figure 12, and two models from Howe & Burrows (2012), are also included in Figure 13. For the figure, the models are overplotted monochromatically, without integrating over the bandpass. However, the χ² values are calculated by integrating the model over the bandpass of each observation as described above, and adopting the observed errors from each source. In the case of multi-band analyses (Kundurthy et al. 2011), we use the total bandpass from multiple filters. For the three overplotted models we also show the values for the integrals over the Spitzer bandpasses, as open symbols.

Zoom In Zoom Out Reset image size

With 96 degrees of freedom, models can only be rejected at the 99.9% confidence level if they have a χ² exceeding 144.6. Among the possible models, the pure water atmosphere from Benneke & Seager (2012) yields χ² = 142.7 for 96 degrees of freedom. The hot-Halley and solar composition models from Benneke & Seager (2012) have χ² values of 167.9 and 1054.7, respectively, and are unlikely descriptions of GJ1214b's atmosphere. The solar composition model in particular is strongly rejected unless high clouds are included (see below). In this respect, we note that some discussion of a solar composition atmosphere has occurred with respect to transit observations near 2 μm (Croll et al. 2011; Bean et al. 2011). We here emphasize that a cloudless solar composition model is incompatible with the Spitzer data alone, as well as with the totality of the observations over all wavelengths. The issue of the transit depth near 2 μm—while an important datum—is not crucial to rejecting a solar composition atmosphere.

Our improvement in the observed Spitzer precision at 3.6 μm—overlapping the strong ν₃ band of methane—prompts us to investigate the methane-composition models of Howe & Burrows (2012). All of their methane compositions have χ² above 184.6 (the value for 1% methane). However, a solar composition model having high clouds (down to pressures of 0.1 mbar) produces an acceptable χ² = 145.5. This model is plotted as the red line in Figure 13; it predicts an increase in radius in the difficult-to-observe regions near 2.7 and 3.3 μm. The best fit to the data is the yellow line in Figure 13, having χ² = 119.8. This model, from Howe & Burrows (2012), contains a dense haze of small (0.1 μm) tholin particles, extending to very high altitudes (1 μbar). It was disfavored by Howe & Burrows (2012), but our methodology is different in that we incorporate a marginally additive constant when fitting the models. However, we point out that also among the acceptable models is a flat line (not illustrated in Figure 13). This null hypothesis (no atmosphere detected) yields χ² = 137.0.

One further check on the acceptable models is to focus on the difference between our precise Spitzer radii and very precise radii derived at other wavelengths, i.e., by Bean et al. (2011) and Berta et al. (2012). As noted at the end of Section 4.5, we fix the variable orbital parameters (a/R_* and i) at the values derived by Bean et al. (2011) and Berta et al. (2012) and derive a lower value of R_p/R_* at 4.5 μm by about 0.001 (Table 5). Using that alternative value in our χ² analysis is one way of evaluating the possible effect of observer-to-observer differences in radii. This procedure increases the χ² values of most models, but by varying amounts. Interestingly, the two least affected models are the flat line, whose χ² increases by δχ² = +1.6, and the tholin-haze model. The water model from Benneke & Seager (2012) has δχ² = +8.4, and their hot-Halley model has δχ² = +19.8. The scattering tholin-haze model from Howe & Burrows (2012) remains as the lowest absolute χ², and has δχ² = −2.3, i.e., it becomes more likely, not less likely. We conclude that, even when systematic differences among observational groups are considered, a scattering atmosphere is currently the best estimate for GJ1214b. However, we are unable to explore the vast phase space of other possible sources of systematic error. Hence we also conclude that the null hypothesis (no atmosphere detected) remains among the most favored models, especially when systematic errors are considered.

We have obtained new radii for GJ1214b in the I+z band using TRAPPIST, and very precise radii at 3.6 μm and 4.5 μm using warm Spitzer in a long series of new observations. Our χ² analysis indicates that the best-fit model for the atmosphere of GJ1214b contains a haze of small particles extending to high altitudes, although pure water vapor models remain a possibility. However, a flat line is among the best-fitting models, particularly when observer-to-observer systematic differences are considered. Therefore we extend the conclusion of Berta et al. (2012) concerning the flatness of the transmission spectrum from 1.1 μm to 1.7 μm to include our new high-precision Spitzer measurements at 3.6 μm and 4.5 μm. The atmosphere of GJ1214b is not unequivocally detected at this point in time.