A SPITZER TRANSMISSION SPECTRUM FOR THE EXOPLANET GJ 436b, EVIDENCE FOR STELLAR VARIABILITY, AND CONSTRAINTS ON DAYSIDE FLUX VARIATIONS

We derive limb-darkening coefficients for the star using a Kurucz ATLAS stellar atmosphere model with T_eff = 3500 K, log (g) = 5.0, and [Fe/H] = 0 (Kurucz 1979, 1994, 2005), where we take the flux-weighted average of the intensity profile in each IRAC band and then fit this profile with four-parameter nonlinear limb-darkening coefficients (Claret 2000). We also derive limb-darkening coefficients for a PHOENIX atmosphere model (Hauschildt et al. 1999) with the same parameters and list both sets of coefficients in Table 3. We trim the maximum stellar radius in the PHOENIX models, which is set to an optical depth of 10⁻⁹, to match the level of the τ = 1 surface in each Spitzer band. We estimate the location of this surface by determining when the intensity relative to that at the center of the star first drops below e⁻¹ and find that the new stellar radius is 0.09%–0.10% smaller than the old τ = 10⁻⁹ value. We find that we can achieve satisfactory four-parameter nonlinear fits to the PHOENIX intensity profiles only when we exclude points where μ < 0.025, whereas the ATLAS models are well described by fits including this region.

As illustrated in Figure 4, the PHOENIX model predicts stronger limb darkening in all bands as compared to the ATLAS model, with the largest differences in the 3.6 μm band. When we compare our best-fit transit parameters using either the ATLAS or PHOENIX limb-darkening coefficients, we find that the best-fit planet–star radius ratios are 0.8σ–1.2σ (0.5%–0.6%) deeper in the 3.6 μm band, 0.06σ–0.07σ (0.04%–0.05%) smaller in the 4.5 μm band, and 0.3σ–0.4σ (0.2%–0.3%) larger in the 8.0 μm band for the PHOENIX models. The best-fit values for the inclination and a/R_⋆ increase by 1.0σ (0.6%) and 0.9σ (0.04%), respectively, for the PHOENIX model fits; this is a product of the stronger limb-darkening profile, as GJ 436b's relatively high impact parameter creates a partial degeneracy between the limb-darkening profile and the other transit parameters.

Zoom In Zoom Out Reset image size

We examine the relative importance of the assumed stellar parameters by comparing two PHOENIX models with effective temperatures of 3400 K and 3600 K. We find that for this 200 K range in effective temperature, the best-fit planet–star radius ratios change by 0.11σ–0.16σ at 3.6 μm , 0.07σ–0.09σ at 4.5 μm,  and 0.10σ–0.12σ at 8.0 μm. The changes in the best-fit values for the inclination and a/R_⋆ were similarly small, 0.04σ and 0.4σ, respectively. We therefore conclude that changes in the stellar effective temperature of less than 200 K are negligible for the purposes of our transit fits. We also compute PHOENIX model intensity profiles for 0.0 < [Fe/H] <+0.3, but we find that the differences between models are much smaller than for our 200 K change in the effective temperature.

As there are currently few observational constraints on limb-darkening profiles for main-sequence stars (e.g., Claret 2008, 2009), and even fewer constraints for M stars in the mid-infrared, we also consider simultaneous transit fits in which we allow quadratic limb-darkening coefficients in each band to vary as free parameters. As a result of the planet's high impact parameter, our observations do not directly constrain the limb-darkened intensity for values of θ ⩽ 50°, corresponding to μ ⩾ 0.64, as the planet does not cross this region on the star. However, we can infer the limb-darkening profile in this region if we assume a simple quadratic limb-darkening law. We require the intensity profile computed from these coefficients to be always less than or equal to one (i.e., no limb brightening) and require the relative intensity at the edge of the star to be greater than or equal to the equivalent K-band limb-darkening from Claret (2000). The dotted lines in Figure 4 show the resulting best-fit limb-darkening profiles in each band; these profiles show less contrast than either model, but the ATLAS models appear to provide the closest match.

This agreement is reflected in the χ² values for the simultaneous transit fits; the total χ² for the best-fit quadratic coefficients is 536,729.25, for the 3500 K ATLAS limb-darkening coefficients it is 536,733.98, and for the [3400, 3500, 3600] K PHOENIX models it is [536,740.69, 536,739.75, 536,738.81], for 536,798 points and either 53 (with fixed limb-darkening) or 59 (with freely varying quadratic limb-darkening coefficients) free parameters. We use the ATLAS limb-darkening coefficients in our subsequent analysis, as they produce a marginally better agreement with the best-fit profiles than the PHOENIX models. Although the χ² value for the best-fit quadratic limb-darkening coefficients is formally smaller than that of either model, this fit also contains six additional degrees of freedom, making the difference negligible.

As an additional test, we also repeat our fits with the limb-darkening coefficients fixed to zero in all bands. This produces planet–star radius ratios that are 1.6σ–2.4σ (1.1%) smaller in the 3.6 μm band, 1.7σ–2.0σ (1.1%) smaller in the 4.5 μm band, and 0.9σ–1.1σ (0.6%–0.7%) smaller in the 8.0 μm band. The best-fit inclination and a/R_⋆ are 2.7σ (1.5%) and 2.0σ (0.1%) smaller, respectively. The χ² value for this fit is 536,738.04, equivalent to the PHOENIX model fits and marginally worse than the ATLAS models or the fitted limb-darkening coefficients. This fit confirms the pattern suggested earlier, namely that stronger limb darkening leads to larger planet–star radius ratios and larger values for the inclination and a/R_⋆. If we consider the constraints imposed by the transit fits, stronger limb darkening means that for a grazing transit the planet must occult a relatively larger fraction of the star in order to produce the same apparent transit depth. This effect will be even larger in visible light, and we conclude that accurate limb-darkening coefficients are essential when calculating the planet–star radius ratio and corresponding transmission spectrum for near-grazing transits.

It is difficult to diagnose the origin of the disagreement between ATLAS and PHOENIX stellar atmosphere models for GJ 436; Kurucz (2005) notes that the ATLAS models should not be used for stellar effective temperatures below 3500 K, as they do not include important low-temperature opacity sources such as TiO and VO. However, these molecules primarily affect the star's visible and near-infrared spectra, and at 3500 K they are still relatively weak (Cushing et al. 2005). Both disk-integrated and intensity spectra for the ATLAS models in this temperature range are featureless longward of 2.4 μm, with the exception of the CO band between 4.3 and 5.0 μm, whereas PHOENIX spectra also show clear molecular band structures, mainly due to H₂O and OH, between 2.5–3.6 μm and 6.5–8.0 μm with corresponding increases in the amount of limb darkening in these bands. The presence of the CO band in both model sets would appear to explain the relatively good agreement in limb-darkening profiles for the 4.5 μm Spitzer bandpass, but we were unable to determine the reason behind the missing mid-IR H₂O absorption features in the ATLAS models, which incorporate the strongest water lines (Kurucz 1999) from the Ames list of Partridge & Schwenke (1997). It is possible that the spherical geometry used in the PHOENIX models (ATLAS models use a plane-parallel geometry) may also affect the resulting limb-darkening profiles (Orosz & Hauschildt 2000; Claret & Hauschildt 2003), but we find that PHOENIX models computed with a planet-parallel geometry show nearly identical limb darkening, with the exception of an exponential decline in the optically thin limb. We conclude that the differing opacities in the 3.6 and 8.0 μm bands appear to be the most likely explanation for the disagreement between the limb-darkening profiles at these wavelengths. In this case, the change in the χ² value indicates that the differences between the two models are not statistically significant for this data set; near-IR grism spectroscopy of transits of GJ 436b, such as those obtained by Pont et al. (2008b), might help to better distinguish between these models.

We calculate uncertainties for our best-fit transit parameters using a Markov Chain Monte Carlo (MCMC) fit (see, for example, Ford 2005; Winn et al. 2007) with a total of 6 × 10⁶ steps, 14 independent chains, and 53 free parameters. Free parameters in the fits include a/R_⋆, i, eight individual estimates of R_P/R_⋆, eight transit times, eight constants, a linear function of time and linear and quadratic terms in x and y for each of the 3.6 and 4.5 μm transits (20 variables total), the amplitude c₂ and decay time c₃ from Equation (2) for the exponential fits to three 8.0 μm transits (6 variables total), and a linear function of time for the other 8.0 μm visit. We assume a constant error for the points in each individual transit light curve, defined as the uncertainty needed to produce a reduced χ² equal to one for the best-fit transit solution.

We initialize each chain at a position determined by randomly perturbing the best-fit parameters from our Levenberg–Marquardt minimization. After running the chain, we search for the point where the χ² value first falls below the median of all the χ² values in the chain (i.e., where the code had first found the optimal fit) and discard all steps up to that point. We calculate the uncertainty in each parameter as the symmetric range about the median containing 68% of the points in the distribution, except for the inclination and a/R_⋆, which we allow to have asymmetric error bars spanning 34% of the points above and below the median, respectively. The distribution of values was very close to symmetric for all other parameters, and there did not appear to be any strong correlations between variables. As a check we also carried out a residual permutation error analysis (Gillon et al. 2007b; Winn et al. 2008), which is sensitive to correlated noise in the light curve, on each individual transit. At the start of each new permutation, we randomly drew values for the inclination and a/R_⋆ from the simultaneous MCMC distribution and then fit for the corresponding best-fit values for the transit time and R_p/R_⋆ in that step. This ensures that our resulting error distributions for individual transit times and R_p/R_⋆ values also take into account the uncertainties in the best-fit values for the inclination and a/R_⋆. In each case where both an MCMC and residual permutation uncertainty are available for a given parameter, we use the higher of the two values. We find that the MCMC fits generally produce larger uncertainties for the 8 μm observations, whereas for 3.6 and 4.5 μm data sets, which have higher levels of correlated noise, the residual permutation uncertainties are typically 50% larger than the MCMC errors.

We fit the secondary eclipses individually using the best-fit values for inclination and a/R_⋆ from our transit fits but allowing the eclipse depths and times to vary freely. Figure 5 shows the final binned data from these fits with the best-fit normalizations for the detector ramp in each channel overplotted, and Figure 6 shows the binned data once these trends are removed, with best-fit eclipse curves overplotted. Figure 7 shows the root-n scaling for the residuals from these fits. Best-fit parameters for individual eclipses are given in Table 4, and the error-weighted average (i.e., weights equal to 1/σ²) of these eclipse depths is reported in Table 5. We find that fixing the time of eclipse to a constant value, defined here as the best-fit orbital phase, does not significantly change our best-fit eclipse depths, nor does it reduce the uncertainties in our measurement of those depths. We calculate the uncertainties on individual eclipses using both an MCMC analysis and a residual permutation error analysis, again taking the higher of the two values as the final uncertainty for each parameter.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We first consider the possibility that the radius of the planet is changing in time, either due to thermal expansion of the atmosphere or the presence of a variable cloud layer at sub-mbar pressures. We require a change in radius of approximately 4% in order to match both of the measured 3.6 μm transit depths; if this change is due to thermal expansion, we can estimate the energy input required using simple scale arguments.

The effective change in the planet's radius due to heating of the atmosphere depends on both the amount of heating and the range in pressures over which this heating takes place. We use the secondary eclipse depths in Section 4.3 to place an upper limit on the allowed change in temperature at the level of the mid-IR photosphere and then calculate the corresponding range in pressure that must be heated by this amount in order to increase the radius of the planet by 4%. If we assume a hydrogen atmosphere with a baseline temperature of 700 K, we find a corresponding scale height of approximately 240 km, where the scale height is defined as $H=\frac{kT}{\mu g}$, T is the temperature of the atmosphere, μ is the mean molecular weight, and g is the surface gravity. We know from the secondary eclipse observations described in Section 4.3 that the temperature of the planet's dayside atmosphere must change by less than 30%, which would correspond to an upper limit of 100 km on corresponding changes in the planet's scale height.

In order to calculate the required energy input to produce the observed change in radius, we must first determine the range of pressures affected by this heating. We model the planet as an interior region with a constant temperature, surrounded by an outer envelope that expands and contracts freely with changing temperature. We set the upper boundary on this region equal to 50 mbar, corresponding to the approximate location of the τ = 1 surface in the mid-infrared. As illustrated in Figure 14, opaque clouds at this pressure suppress but do not entirely remove absorption features in the planet's transmission spectrum at these wavelengths, making this a reasonable estimate for the location of the τ = 1 surface. We assume that when the planet is heated the scale height changes by 100 km, which requires the lower boundary of the heated region to be located at a pressure of approximately 1 bar in order to produce a 1% expansion in radius. If we then calculate the change in the planet's gravitational energy corresponding to this expansion, we find that an energy input of approximately 10²⁶ J is required. Repeating this calculation for a 4% increase in radius, we find a lower boundary at 8000 bars and a corresponding energy input of 10³⁰ J. The insolation received by the planet is 10²⁰ W, which gives an energy budget of 10²⁵ J per orbit. When we examine Figure 10 we find that the observed change in radius occurs primarily between the third and fourth visits (UT 2009 January 25–28). This would require an energy input as much as 10⁵ times higher than the total insolation over this epoch, which is clearly unphysical.

One alternative explanation for the observed change in radius would be to invoke the presence of intermittent, high-altitude clouds. Such clouds could produce a change in the apparent radius of the planet across multiple bands without requiring any actual heating or cooling of the atmosphere. In this picture, smaller radii for the planet would correspond to the cloud-free state, while larger radii would require the presence of an additional cloud layer. A change of 4% in apparent radius would require the clouds to form at a pressure approximately 100 times lower than the location of the nominal cloud-free radius. In Section 4.2, we find that the average pressure of the τ = 1 surface for the nominal methane-poor (green) model between 3 and 10 μm is 40 mbar, indicating that the clouds would have to extend to 0.4 mbar to explain the largest measured 3.6 μm radius for the planet. This conclusion is reasonably independent of our assumed composition, as the average τ = 1 surface for the methane-rich (blue) model is located at 30 mbar. Gravitational settling would presumably pose a challenge for cloud layers at sub-mbar levels, but vigorous updrafting of condensate particles might compensate for this effect. The broadband nature of the data presented here makes it difficult to directly test this hypothesis; we therefore recommend the acquisition of high signal-to-noise, near-infrared grism spectroscopy over multiple transits in order to resolve this issue. A 0.5 mbar cloud layer would lead to a near-featureless transmission spectrum whereas a lower cloud layer would still exhibit many of the same absorption features as a cloud-free atmosphere. Such a data set would also allow us to test the theory, outlined in Section 4.1.3, that the observed transit depth variations are due to the occultation of regions of non-uniform brightness on the surface of the star, as these regions should also produce a wavelength-dependent effect.

It is possible that poorly corrected instrument effects, such as the intrapixel sensitivity variations at 3.6 and 4.5 μm, or the detector ramp at 8.0 μm, might lead to variations in the measured transit depth. Because there is complete overlap between the positions spanned by the star in the in-eclipse and out-of-eclipse data for all 3.6 and 4.5 μm visits, fits that inadequately describe the pixel response as a function of position should fail equally for both sections of the light curve. The UT 2009 January 30 transit serves as an example of imperfectly removed detector effects, as the residuals display a sawtooth signal with a shape and timescale similar to the original intrapixel sensitivity variations (see Section 2.1 for a more detailed discussion of this light curve). Conversely, it is much more difficult to explain the 3.6 μm transit on UT 2009 January 28 with this scenario, as there appears to be a large dip in residuals during ingress, but when the star spans the same pixels in the out-of-eclipse data we see no comparable deviations.

In this section, we consider an alternate decorrelation function that better accounts for small-scale variations in the intrapixel sensitivity function as discussed in Ballard et al. (2010b). Following the discussion in Ballard et al., we describe the intrapixel sensitivity variations using a position-weighted average of the time series after the best-fit transit function and linear function of time from the position fits described above has been divided out. Unlike Equation (1), this formalism does not assume a functional form for the intrapixel sensitivity variations and therefore should in principle produce an unbiased correction for these variations. We calculate the weighting function as follows:

$$\begin{eqnarray} &&W(x_i,y_i)=\frac{\sum \nolimits _{i\ne j}{\exp \left(-\frac{(x_j-x_i)^2}{2\sigma _x^2}\right)\exp \left(-\frac{(y_j-y_i)^2}{2\sigma _y^2}\right)f_j}}{\sum \nolimits _{i\ne j}{\exp \left(-\frac{(x_j-x_i)^2}{2\sigma _x^2}\right)\exp \left(-\frac{(y_j-y_i)^2}{2\sigma _y^2}\right)}},\quad \end{eqnarray} \tag{ 4 }$$

where x_i and y_i are the x and y positions of the ith frame, x_j and y_j are the x and y positions for the rest of the time series. We optimize our choice of σ_x and σ_y to produce the smallest possible scatter in the final time series when we fix the transit light curve to the best-fit solutions listed in Table 2. We find that the preferred values range between 0.0053–0.0120 pixels in σ_x and 0.0024–0.0045 pixels in σ_y for the four 3.6 and 4.5 μm transits examined here. For ease of computation, we bin our time series in intervals corresponding to one point per original set of 64 images (in some instances there are less than 64 images in a given bin after removing outliers) and iteratively calculate the weighting function and the linear function of time plus transit fits until we converge to a consistent solution.

Once we have a final solution we calculate the weighting function for the unbinned data and carry out a final fit for the transit function to determine our best-fit transit depth. In this case, we fix the inclination and a/R_⋆ to their best-fit values from the simultaneous fits to all transits described in Section 2.3, which allows us to fit each transit individually using the weighting function while still preserving the constraints imposed in a simultaneous fit. We find that in all cases we obtain transit depths and times that are consistent with the values from our fits using Equation (1), with a standard deviation that is comparable or slightly worse than that achieved with our polynomial fits.

We also carried out a second set of fits in which we derived our corrections for the intrapixel sensitivity variations using only the out-of-transit data and found that our best-fit planet–star radius ratios changed by less than 0.4σ in all cases. Because the star samples the same regions of the pixel in both the in-transit and out-of-transit data, it is possible to obtain an equivalently good correction for the intrapixel sensitivity variations using only the out-of-transit points. Conversely, this means that poor corrections for this effect should produce equally large deviations in both the in-transit and out-of-transit regions of the light curve. As we will discuss in the following section, we find that the residuals for the deepest transits in these two bands have a significantly higher rms in transit than out of transit. This behavior is inconsistent with our expectations for poorly corrected instrument effects, and we therefore conclude that it is unlikely that these effects are responsible for the discrepant transit depths measured at 3.6 and 4.5 μm.

At 8.0 μm we fit the data with a single or double exponential function to describe the smoothly varying detector ramp. In Agol et al. (2010), we conclude that this functional form avoids correlations between the slope of the ramp and the measured transit or eclipse depth; however, we check this assertion using our 8 μm data as well. For our 8 μm transit fits, we find that the exponential term has a coefficient of [0.00156, 0.00000, 0.00288, 0.00299], corresponding to planet–star radius ratios of [0.08234, 0.08162, 0.08138, 0.08336], where we have set the amplitude of the exponential term to zero for the transit occurring in the middle of our 70 hr phase curve observation. For the 11 secondary eclipse observations, we find coefficients of [0.00645, 0.00627, 0.00194, 0.00433, 0.00534, 0.00140, 0.00000, 0.00359, 0.00321, 0.00353, 0.00262], corresponding to eclipse depths of [0.0552, 0.0507, 0.0395, 0.0495, 0.0367, 0.0523, 0.0421, 0.0386, 0.0491, 0.0397, 0.0441], respectively, where we have set the exponential coefficient to zero for the secondary eclipse at the end of the phase curve observation. We find no evidence for any correlation between the slope of the exponential function and the measured transit or secondary eclipse depths. As an additional check we also confirm that there is no correlation between these depths and either the measured sky background or the total stellar flux given in Table 1.

The presence of spots or faculae on the visible face of the star can have two distinct effects on the measured light curve for a transiting planet. Non-occulted spots on the visible face of the star reduce the star's total flux, increasing the measured transit depth, while spots occulted by the planet cause a small positive deviation in the light curve with a timescale proportional to the physical size of the occulted spot (e.g., Rabus et al. 2009) and occulted faculae would have the opposite effect. The early K dwarf HD 189733 (T_eff = 5100 K) is perhaps the best-studied example of an active star with a transiting hot Jupiter (e.g., Bakos et al. 2006; Pont et al. 2008a; Désert et al. 2011a), but the late G dwarf CoRoT-2 (T_eff = 5600 K) also exhibits a high level of spot activity that may have resulted in early overestimates of its planet's inflated radius (Guillot & Havel 2011). This problem is likely to be even more common for M dwarfs, and in fact several instances of occulted spots were reported in transit light curves for the super-Earth GJ 1214b, which orbits a 3000 K primary (Carter et al. 2011; Berta et al. 2011; Kundurthy et al. 2011).

Although it is important to correct for these effects in any transit fit, it is particularly crucial when comparing non-simultaneous, multi-wavelength transit observations such as the ones described in this paper, which require a relative precision of better than a part in 10⁻⁴ in the measured transit depths. We evaluate the likely impact of GJ 436's activity on the measured transit depths using several complementary approaches. First, we estimate the average activity level on GJ 436 by measuring the amount of emission in the cores of the Ca ii H&K lines; in Knutson et al. (2010), we determined that GJ 436 had an average S_HK of 0.620. Isaacson & Fischer (2010) found that other stars in the California Planet Search database with similar B − V colors have S_HK values ranging between 0.5 and 2.0, indicating that GJ 436b is relatively quiet for its spectral type. Demory et al. (2007) report that this star's rotation period is greater than 40 days, consistent with upper limits on vsin i of 1 km s⁻¹ (Jenkins et al. 2009), also suggesting that it is likely to be relatively old and correspondingly quiet. The upper limit of 3 km s⁻¹ on vsin i from spectroscopy (Butler et al. 2004) is also consistent with an inclined or pole-on viewing geometry, although it is not required as long as the star's rotation period is longer than 7 days.

Rather than relying on these indirect measures of activity, we can also directly measure the amplitude of the star's rotation-modulated flux variations using visible-light ground-based observations. We obtained observations of GJ 436 in Strömgren b and y filters over a span of approximately six months surrounding our 2009 Spitzer transit and secondary eclipse observations from an ongoing monitoring program carried out with the T12 0.8 m Automatic Photoelectric Telescope (APT) at Fairborn Observatory in southern Arizona (Henry 1999; Eaton et al. 2003; Henry & Winn 2008). In these observations, the telescope nodded between GJ 436 and three comparison stars of comparable or greater brightness, which were then used to correct for the effects of variable seeing and airmass. We find that during the period between UT 2009 January 9 and February 4, when a majority of our transit data were obtained, the star varied in flux by less than a few mmag in visible light (Figure 11). We carry out a similar check for variability in the infrared using the fifteen 8 μm flux estimates listed in Table 1, which we find have a standard deviation of 0.07%. Both of these measurements indicate that the star is very nearly constant in flux in both visible and infrared light, and we can therefore rule out non-occulted spots as the cause of the observed transit depth variations.

Zoom In Zoom Out Reset image size

We also use these same data to search for periodicities corresponding to GJ 436b's rotation period. If we fit the combined b and y band fluxes with a sine function plus a quadratic function of time as shown in Figure 11, we find a best-fit period of 56.5 days. We calculate a Lomb–Scargle periodogram (Lomb 1976; Scargle 1982) for these data and find that this period has a false alarm probability of only 2%, which we determine using a bootstrap Monte Carlo analysis. We find a nearly identical best-fit period of 56.6 days in the S_HK values measured with Keck HIRES during this epoch (Isaacson & Fischer 2010), but with a much higher false alarm probability of approximately 20%. We also examine the correlation between the measured b fluxes and S_HK values over the six-year period in which both were available (Figure 12) and find that these parameters are negatively correlated. Taken together, these data indicate that the small observed variations in GJ 436's visible-light fluxes are likely connected with the presence of regions of increased magnetic activity on the visible face of the star.

Zoom In Zoom Out Reset image size

Although such low-amplitude flux variations generally indicate that a star has relatively few spots, there are two important exceptions. First, if the spots are uniformly distributed in longitude, it is theoretically possible to have a star with significant spot coverage and an effectively constant flux. It would not be surprising if the occurrence rate and distribution of spots was different for M stars than for G or K stars, but in GJ 436b's case the lack of any flux variations larger than a few mmag would seem to place a strong limit on the allowed spot distributions. We can quantify this limit if we assume that the deviation of approximately 0.08% in the first part of the 3.6 μm transit light curve from UT 2009 January 28 shown in Figure 2 is due to the occultation of a bright region on the star. This region must have a surface intensity that is 12% brighter than the rest of the star in order to produce the observed deviation. If we compare PHOENIX models with varying effective temperatures integrated over this band, we find that the star's temperature must increase by approximately 200 K in the affected region in order to match this surface intensity. We know that the total rotational modulation in the star's visible-light flux must remain below 0.1%, and we estimate that an increase of 12% in the 3.6 μm surface intensity should produce an increase of approximately 65% in the Strömgren (b + y)/2 band. In this case, the fractional area covered by active regions on the star must vary by less than 0.15% from the most active to the least active hemisphere. Of course, it is possible that the stellar atmosphere models do not provide an accurate match for the spectra of these active regions; if we instead use the measured 3.6 μm flux contrast of 12%, we find a more conservative limit of 1% on variations in the area affected by stellar activity.

A second, more plausible scenario involves tilting the rotation axis of star so that we are viewing it closer to pole-on, which would effectively suppress the amplitude of rotational flux variations regardless of spot coverage. If we assume that the star's spin axis is randomly oriented with respect to our line of sight, the probability that it will fall within 45° of a pole-on view is 30%. In this scenario, the star could be highly spotted, allowing for frequent occultations of spots by the planet, while still displaying a small rotational flux modulation. This scenario would require the planet's orbit to be misaligned with respect to the star's rotation axis, but such misalignments are commonly seen in other transiting planet systems (Winn et al. 2010a). Although the Rossiter–McLaughlin effect has never been successfully measured for GJ 436b, Winn et al. (2010b) find that the Neptune-mass planet HAT-P-11b, which is perhaps the best analogue to the GJ 436 system, has a sky-projected obliquity of $103\mbox{$^\circ $}^{+26\mbox{$^\circ $}}_{-10\mbox{$^\circ $}}$ indicating that this system is significantly misaligned. If most close-in planets start out misaligned and are then gradually brought into alignment through tidal interactions with their host star as proposed by Winn et al. (2010a), the fact that HAT-P-11b still maintains both a non-zero orbital eccentricity and a significant misalignment would seem to suggest that the same could also be true for GJ 436b.

If we proceed with the hypothesis that GJ 436 is both spotty and tilted with respect to our line of sight, we can then search for evidence of occulted spots in the light curves with discrepant transit depths. We first compare the relative standard deviations of the in-transit ($\sigma _{\tt in}$) and out-of-transit ($\sigma _{\tt in}$) residuals plotted in Figure 2:

$$\begin{eqnarray} &&\sigma _{\tt rel} = \frac{\sigma _{\tt in} - \sigma _{\tt out}}{\sigma _{\tt out}}. \end{eqnarray} \tag{ 5 }$$

We list the measured values of $\sigma _{\tt rel}$ for all eight transit observations in Table 7. Both the 3.6 μm transit on UT 2009 January 28 and the 4.5 μm transit on UT 2009 January 30 appear to have inflated values of $\sigma _{\tt rel}$, as would be expected if the planet occulted active regions on the star during these visits. We can quantify the statistical significance of the measured $\sigma _{\tt rel}$ values if we assume that both the in-transit and out-of-transit points are drawn from the same underlying Gaussian distribution and then ask how many times in a sample of 100,000 random trials we measure a value of $\sigma _{\tt rel}$ greater than or equal to the value calculated directly from our observations. In each trial, we generate two synthetic data sets, each with the appropriate length corresponding to either the in-transit or out-of-transit measurements, and then calculate the standard deviation of each distribution and the corresponding value of $\sigma _{\tt rel}$. In the 3.6 μm transit observation on UT 2009 January 9, there are 81,848 out-of-transit flux measurements and 25,482 in-transit flux measurements, and we find that over 100,000 trials, we obtain a value of $\sigma _{\tt rel}$ greater than or equal to the measured value of 0.2% approximately 36% of the time. Repeating the same calculation for the 3.6 μm transit observed on UT 2009 January 28, which has 82,238 out-of-transit points and 25,530 in-transit points, we obtain $\sigma _{\tt rel}$ greater than or equal to the measured value of 1.4% only 0.23% of the time. We list the corresponding probabilities for all eight transits in Table 7.

Table 7. A Comparison of the In-transit versus Out-of-transit Standard Deviations

UT Date	λ (μm)	$N_{\tt in}^{\rm a}$	$N_{\tt out}^{\rm a}$	$\sigma _{\tt rel}$	$P(\sigma _{\tt rel})^{\rm b}$
Unbinned data
UT 2007 Jun 29	8.0	7,924	17,956	−1.3%	0.92
UT 2008 Jul 14	8.0	7,895	25,012	+1.4%	0.059
UT 2009 Jan 9	3.6	25,482	81,848	+0.2%	0.36
UT 2009 Jan 17	4.5	25,954	82,212	−0.3%	0.70
UT 2009 Jan 25	8.0	7,910	25,334	+0.1%	0.44
UT 2009 Jan 28	3.6	25,536	82,238	+1.4%	0.0023
UT 2009 Jan 30	4.5	25,955	62,334	+1.1%	0.018
UT 2009 Feb 2	8.0	7,890	25,318	+0.1%	0.45
Binned data
UT 2007 Jun 29	8.0	126	287	−13.6%	0.97
UT 2008 Jul 14	8.0	126	399	−4.9%	0.75
UT 2009 Jan 9	3.6	411	1311	+3.3%	0.21
UT 2009 Jan 17	4.5	412	1310	−3.4%	0.80
UT 2009 Jan 25	8.0	126	403	+9.7%	0.093
UT 2009 Jan 28	3.6	411	1311	+37.5%	1 × 10−6
UT 2009 Jan 30	4.5	412	989	−1.4%	0.63
UT 2009 Feb 2	8.0	126	403	+2.9%	0.34

Notes. ^aNumber of in-transit and out-of-transit points. ^bProbability that the standard deviation of the in-transit data would be greater than the standard deviation of the out-of-transit data by an amount $\sigma _{\tt rel}$ if both data sets are drawn from the same underlying Gaussian distribution.

Download table as:  ASCIITypeset image

We also repeat this same test with data that have been binned in sets of 64 images, corresponding to 10 s bins at 3.6 and 4.5 μm and 30 s bins at 8 μm. This allows us to evaluate the relative contribution that correlated noise makes to the in-transit and out-of-transit variances, as the photon noise should be reduced by a factor of eight in these bins (also see Figure 3). In this case, we carry out 1,000,000 random trials for each visit, as each simulated data set is much smaller and the computations are correspondingly fast. We find that for the binned January 9 light curve there are 1311 points out of eclipse and 411 points in eclipse. In this case, $\sigma _{\tt rel}$ is 3.3%, and we obtain values greater than or equal to this number in 21% of our random trials. Repeating this calculation for the UT 2009 January 28 visit, we find that the measured value of $\sigma _{\tt rel}$ is 37% (i.e., a standard deviation that is 37% higher in eclipse than it is out of eclipse), with 1311 points out of eclipse and 411 points in eclipse. In our simulations, assuming a single Gaussian probability distribution for both segments, this level of disagreement occurred only once in 10⁶ trials. We find that in all other visits, including the 4.5 μm transit observed on UT 2009 January 30, the binned data in and out of eclipse are consistent with a single distribution.

One consequence of a misalignment between the star's rotation axis and the planet's orbit is that the planet will not necessarily occult the same spot on successive transits, as would be expected for a well-aligned system; we therefore consider each transit individually. Our analysis above indicates that the 3.6 μm transit on UT 2009 January 28 displays a statistically significant increase in the standard deviation of the in-transit data that is dominated by contributions from correlated noise on timescales greater than 30 s, as would be expected if the planet occulted an active region on the surface of the star. Although the 4.5 μm transit from UT 2009 January 30 does not appear to display a similar increase, our imperfect correction for the intrapixel sensitivity variations in this visit means that we are less sensitive to variations in $\sigma _{\tt rel}$. We argue that even if the star's rotation axis and the planet's orbit are misaligned, it is still likely that the planet would occult the same active region during both the UT 2009 January 28 and January 30 visits, as the interval between these visits is much shorter than the star's approximately 50 day rotation period. As we discuss later in this section, the fact that both visits display increased transit depths and shifted transit times provides additional support for this hypothesis.

We also consider the possibility that the increased scatter in the in-transit residuals might be due to a change in the transit parameters, including the planet's radius, orbital inclination, transit time, or a/R_⋆, from one visit to the next. We test this hypothesis by taking the difference of the first and second visits in each bandpass from 2009 and comparing the shape of the residual light curve to the differences we would expect due to changes in these parameters, which should be distinct from the deviations created by occulted star spots (Figure 13). Because we are directly differencing the two light curves, our results are independent of any assumptions about the shape of the transit light curve or the stellar limb darkening. We inspect the deviations in the residuals plotted in Figure 13 and conclude that they do not appear to be well matched by changes in the best-fit transit parameters, leaving occultations of active regions on the surface of the star as the most likely hypothesis.

Zoom In Zoom Out Reset image size

If the planet occults a spot it can also cause a shift in the best-fit transit times, particularly when the spot is near the edge of the star and is occulted during ingress or egress. Indeed, we see that the UT 2009 January 28 3.6 μm appears to occur 31.4 ± 9.5 s early, while the 4.5 μm January 30 visit occurs 34.4 ± 9.4 s late (see Figure 8) in the fits where we fix a/R_⋆ and i to a single common value. As a test we repeated our fit to the 3.6 μm transit excluding the first 1/3 of the transit light curve and found that the best-fit transit time shifted forward by approximately 30 s. We would also expect that transits observed in visible light, where the contrast between the spots and the star is more pronounced, would show proportionally larger timing deviations when the planet crosses a spot. As noted in Section 3.1, the scatter in the measured visible-light transit times is inconsistent with a constant period, and the amplitude of the visible-light deviations is on average larger than the deviations in the infrared. We should also see this same wavelength dependence in the measured transit depths in Figure 10, and indeed we find that the 3.6 μm transit depth changes by 7.8%, the 4.5 μm transit depth changes by 5.3%, and the 8.0 μm transit depth changes by 4.9% during the period between UT 2009 January 9 and February 2. Lastly, we can examine the visible-light flux measurements for GJ 436 in Figure 11 and see that these two transits were obtained near a minimum in the star's flux, consistent with a relative increase in the fractional spot coverage as compared to earlier epochs. The measured values for S_HK, a common activity indicator, appear to be anti-correlated with the observed flux variations and reach a local maximum near this point.

Tidal dissipation is expected to have driven GJ 436b into a pseudosynchronous rotation state in which the planet's spin frequency is nearly commensurate with the planet's instantaneous orbital frequency at periastron. There are several competing theories of the pseudosynchronization process (see, e.g., Ivanov & Papaloizou 2007). We adopt the expression given by Hut (1981):

$$\begin{eqnarray} &&{\Omega _{\rm spin}\over {\Omega _{\rm orbit}}}={ 1+{15\over {2}}e^{2}+{45\over {8}}e^{4}+{5\over {16}}e^{6} \over {\big(1+3e^2+{3\over {8}}e^{4}\big)(1-e^{2})^{3/2}}}{.} \end{eqnarray} \tag{ 6 }$$

For GJ 436b, this relation gives P_spin = 2.32 days, which yields a 19 day synodic period for the star as viewed from a fixed longitude on the planet. GJ 436b also experiences an 83% increase in incident flux during the 1.3 day interval between apastron and periastron.

We have computed simple hydrodynamical models to assess whether the asynchronous rotation and time-varying insolation are likely to generate atmospheric flows that are sufficiently chaotic to produce observable orbit-to-orbit variability in the secondary eclipse depths. Our two-dimensional hydrodynamical model contains three free parameters. The first, p_8 μm, is the atmospheric pressure at the 8 μm photosphere; the second, X, corresponds to the fraction of the incoming optical flux that is absorbed at or above the 8 μm photosphere; and the third, p_b, corresponds to the pressure at the base of our modeled layer. We adopt parameter values of p_8 μm = 100 mbar, p_b = 4.0 bar, and X = 1.0 for these models, which put our model's light curve in good agreement with GJ 436b's average 8 μm secondary eclipse depth. The full details of the computational scheme are the same as those adopted in Langton & Laughlin (2008), with updates as described in Laughlin et al. (2009). A model photometric light curve is then obtained by integrating at each time step over the planetary hemisphere visible from Earth, where we assume that each patch of the planet radiates with a blackbody spectrum corresponding to the local temperature.

The model is run for a large number of orbits, and a quasi-steady state surface flow emerges. The temperature structure of this flow as seen from an observer in the direction of Earth at five equally spaced intervals in the orbit is shown in Figure 17, and the model light curve over these five orbits is shown in Figure 18. Over the course of a single orbit, the 8 μm planet-to-star flux ratio varies nearly sinusoidally from ΔF/F = 0.033% to 0.043%. The model's flux at secondary eclipse agrees well with the observed value and varies by only 0.5% peak to peak from one orbit to the next. We note that more sophisticated three-dimensional general circulation models for GJ 436b from Lewis et al. (2010) also predict very low (1.3%–1.5%) levels of variability in the 8 μm band for a range of atmospheric metallicities.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Although these models indicate that GJ 436b's modest orbital eccentricity is likely not sufficient to induce significant variability, they also do not include many processes such as clouds, photochemistry, and small-scale turbulence that are known to contribute to temporal variability in planetary atmospheres. We therefore place empirical limits on GJ 436b's dayside variability using the eleven 8 μm secondary eclipse observations. We assume that the intrinsic dayside fluxes are drawn from either a Gaussian distribution with a standard deviation δ or from a boxcar distribution with a width equal to 2δ. In both cases, we set the mean of the distribution equal to the error-weighted mean of the measured secondary eclipse depths given in Table 5. We then conduct 10,000 random trials, where we draw 11 measurements from each distribution and calculate the reduced χ² of these values as compared to the measured secondary eclipse depths in Table 3. We then determine the fraction of the 10,000 random trials in which the reduced χ² is less than or equal to one, which should correspond to the probability that the underlying distribution is consistent with the measured eclipse depths. We repeat this calculation for a range of values for δ and plot the resulting probability distribution as a function of δ for both boxcar and Gaussian distributions. We find that for a boxcar distribution we can place [1σ, 2σ, 3σ] limits on the intrinsic variability of [29%, 42%, 58%], and for a Gaussian distribution our corresponding upper limits are [17%, 27%, 42%]. These limits are consistent with the predictions from general circulation models for this planet, but they are not low enough to provide meaningful constraints on these models.