THE ATMOSPHERES OF THE HOT-JUPITERS KEPLER-5b AND KEPLER-6b OBSERVED DURING OCCULTATIONS WITH WARM-SPITZER AND KEPLER

As Spitzer exhausted its cryogen of liquid coolant on 2009 May 15, only the first two channels, at 3.6 (channel 1) and 4.5 μm (channel 2), of the Infrared Array Camera (IRAC; Fazio et al. 2004) are available in the post-cryogenic mission. Kepler-5 and Kepler-6 were both visited twice at each available IRAC bandpass, leading to a total of eight eclipses being secured as part of the 800 hr allocated to the program PID 60028 (PI: D. Charbonneau). We present here the first observations from this program which contribute to characterizing and vetting Kepler candidates. Each visit is secured in full frame mode (256 × 256) with an exposure time of 10.4 s, at 12 s cadence, leading to 2700 frames per Kepler-5 (K = 11.77) 10 hr observations and 2150 frames per Kepler-6 (K = 11.71) 8 hr observations. The full set of Warm-Spitzer observations is presented in Table 2.

We use the basic calibrated data frames produced by the standard IRAC calibration pipeline for the photometric extraction. These files are corrected for dark current, flat-fielding, and detector non-linearity and converted into flux units.

The first step of the photometric extraction consists of determining the centroid position of the stellar point-spread function (PSF) using DAOPHOT-type Photometry Procedures, CNTRD, from the IDL Astronomy Library.¹³ We use the APER routine to perform an aperture photometry with a circular aperture of variable radius, using radii of 1.5 to 6 pixels, in 0.5 steps. Finally, the propagated uncertainties are derived as a function of the aperture radius, and we adopt the aperture which provides the smallest errors. We find that the eclipse depths vary by less than 1.5σ with the aperture radius for all the light curves analyzed in this project. In the case of Kepler-5, the optimal apertures are around 5 and 2.5 pixels at 3.6 and 4.5 μm, respectively. Kepler-6 has a companion at an angular separation of 41 and which is 3.8 mag fainter (Dunham et al. 2010). This star is located at 3.5 IRAC pixels from our target of interest. Therefore, we fix our circular aperture at a radius of 2.5 pixels for all the visits of this star in order to minimize the contribution of the companion.

The background level for each frame is determined by two methods. We use APER to measure the median value of the pixels inside an annulus centered on the star with inner and outer radii of 9 and 16 pixels, respectively. We also estimate the background by fitting a Gaussian to the central region of the histogram of counts from the full array. The center of the Gaussian fit is adopted as the residual background intensity. These two methods provide the same background level; therefore, we use the background values determined from an annulus.

The contribution of the background to the total flux from the stars Kepler-5 and Kepler-6 is low in both IRAC bandpasses, from 0.2% to 0.55% depending of the photometric aperture size. Therefore, photometric errors are not dominated by fluctuations in the background. We note that the background is negative for all channel 1 observations and for one observation secured in channel 2. The 10.4 s exposures yield a typical signal-to-noise ratio of 200 and 150 per individual observation at 3.6 and 4.5 μm, respectively.

After producing the photometric time series, we use a sliding median filter that compares the 20 preceding and 20 following photometric measurements and centroid positioning to identify and trim outliers greater than 5σ. This process removes measurements affected by transient hot pixels and inaccurate centroid determination. In this way, we discarded approximately 2% photometric points from all the observations. We also discarded the first half hour of observations, corresponding to around a hundred frames. These frames exhibit an anomalously large pointing drift, most likely due to settling of the telescope at the new position. The final number of photometric measurements used for each observation is given in Table 2.

To measure the occultation depths and their uncertainties we model the light curves with four parameters: the occultation depth d, the orbital semimajor axis to stellar radius ratio (system scale) a/R_⋆, the impact parameter b, and the time of mid transit T_c. We use the IDL transit routine OCCULTSMALL, developed by Mandel & Agol (2002), to model the light curve. Kepler-5b and Kepler-6b have both well-defined transit and stellar parameters constrained by the transits from the Kepler light curves (Koch et al. 2010a; Dunham et al. 2010). Thus, we adopt these parameters, in particular the inclination and the scale of the system, by fixing a/R_⋆ and b to their nominal values in our eclipsing model. We also fixed the mid-eclipse times to the values we derived from the Kepler light curves (see Section 3 below). Only the depth of the occultations are allowed to vary when fitting the Warm-Spitzer observations.

The Warm-Spitzer/IRAC photometry is known to be systematically affected by the so-called pixel-phase effect, which is due to the combination of pointing jitter and intra-pixel sensitivity (see, e.g., Charbonneau et al. 2005, 2008; Reach et al. 2005; Désert et al. 2009, 2011). This effect corresponds to oscillations in the measured raw light curve with an approximate period of 70 minutes and an amplitude of 2% peak-to-peak. This artifact has to be corrected to properly extract the eclipse depths and errors. To correct the raw light curve for this intra-pixel sensitivity effect, we use the centroid position of the target on the detector and its variations as function of time. We use a quadratic function of the position and a linear function of time F_corr = F[K₁(x − x₀) + K₂(x − x₀)² + K₃(y − y₀) + K₄(y − y₀)² + K₅ + K₆ × t], where F and F_corr are the fluxes of the star before and after the pixel-phase effect correction, and (x − x₀) and (y − y₀) define the position in pixels of the source centroid on the detector with respect to the pixel pointing position, located at [x₀, y₀], and the constants K_i are the free parameters to adjust.

We use the MPFIT package¹⁴ to perform a Levenberg–Marquardt least-squares fit of the transit model to observations. The best-fit model is computed over the whole parameter space (d, K_i). The baseline function described above is combined with the transit light curve function so that the fit is constrained by eight parameters (two for the transit model, two for the linear baseline, and four for the pixel phase effect).

We use four methods to determine the centroid position of a stellar PSF and test the its impact on the final results. In the first method, we fit Gaussians to the marginal x and y sums using GCNTRD. The second method we apply consists of computing the centroid of the stars using a derivative search with CNTRD also part of the standard IDL Astronomy Library. As for the third method, we calculate the center-of-light of the star within a circular aperture with a radius of 3.0 pixels to approximate the center of the star. Finally, the fourth method uses a symmetric two-dimensional Gaussian fit with a fixed width to a 3 × 3 pixel sub-array to approximate the center of the star, as suggested by Agol et al. (2010).

All these methods allow us to measure the centroid position of a stellar image with a different level of accuracy and precision as previously noted by Agol et al. (2010). For all the methods tested, we find that the centroid position varies by less than 20% of a pixel during a complete observation and that can be determined to a precision of a hundredth of a pixel. We measure the eclipse depths and associated errors for all the four methods following the algorithm we describe below. Although each method provides slightly different results for the centroid determination, it yields an eclipse depth consistent to within 1σ and similar errors. In this paper, we use the CNTRD program as it produces the smallest reduced χ² .

As shown by Pont et al. (2006), the existence of low-frequency correlated noise (red noise) between different exposures must be considered to obtain a realistic estimation of the uncertainties. To obtain an estimate of the systematic errors in our observations we use the permutation of the residuals method, known as "prayer-bead" (Moutou et al. 2004; Gillon et al. 2006, 2007) and derived the covariance from the residuals of the light curve. In this method, the residuals of the initial fit are shifted systematically and sequentially by one frame, and then added to the transit light curve model before fitting again. The error on each photometric point is the same and is set to the rms of the residuals of the first best fit obtained. We find that the final rms values are 5% to 20% larger than the predicted photon noise depending on the target and the bandpass.

We produce as many shifts and fits of transit light curves as the number of photometric measurements to determine the statistical and systematical errors for the adjusted parameters. We set our uncertainties equal to the range of values containing 68% of the points in the distribution in a symmetric range about the median for a given parameter. We check that these values are close to the standard deviation of each nearly Gaussian distribution. We fit for the eclipse depth, fixing the mid-eclipses time to the values we derived from the Kepler observations as describe in the following section. The transit light curves are presented in Figure 1 and the best fits of the occultation depths and their error bars are listed in Table 2.

Zoom In Zoom Out Reset image size

The current upper limit on the orbital eccentricity e from radial velocity measurements is consistent with zero for both targets (Koch et al. 2010a; Dunham et al. 2010). An independent analysis of the two first Kepler quarters (Q0 and Q1) have already been applied (Kipping & Bakos 2011). It uses the transits, eclipses, and radial velocity measurements to confirm that no eccentricity is found for Kepler-5 and to report a marginal eccentricity at 2σ level for Kepler-6.

We measure the mid-occultation timing offset from both Kepler and Warm-Spitzer observations. The determination of the timing of the secondary eclipse constrains the planet's orbital eccentricity. A non-zero value of e could produce a measurable shift in mid-eclipse times. The Warm-Spitzer observations are affected by systematics (correlated noise) which prevent us from deriving timing measurements with a better precision than the precision we derive from Kepler measurements. Indeed, we measure the mid-occultation timing offset from Kepler with uncertainties of nearly 8 minutes for both targets. Our estimate for the best-fit timing offset translates to a constraint on e and the argument of pericenter ω. The timing is used to constrain the ecos (ω) at the 3σ level. For Kepler-5b we find ecos (ω) = −0.025 ± 0.005, and |ecos (ω)| < 0.06 to 3σ. For Kepler-6b we find ecos (ω) = −0.01 ± 0.005, and |ecos (ω)| < 0.035 to 3σ. These upper limits imply that the orbits of these objects are nearly circular unless the line of sight is aligned with the planet's major axes, i.e., the argument of periapse ω is close to 90° or 270°. Assuming that both planets are on a circular orbit and accounting for the time taken by the light to cross the entire orbits, we would expect that the mid-secondary eclipses are delayed by 50 s and 45 s for Kepler-5 and Kepler-6, respectively, which translate into orbital phase of 0.50016 for both targets.

The occultation depths measured in each bandpass are combined and are turned into an emergent spectrum for each planet. The observed flux of the planet in each bandpass corresponds to the sum of the reflected light, the thermal emission of the incident stellar flux, and the interior flux from the planet itself (e.g., internal heat or emission due to tidal forces). As a first step, we assume that the thermal emission from the planet itself is negligible. We use the chromatic information to distinguish among the different components of the observed planetary flux.

If we were to explain the occultation detections in the Kepler bandpass as reflected light for both objects, it would imply that we are indeed measuring their geometric albedos. One can estimate the reflected light by the planet as F_p/F_⋆ = A_g(R_p/a)², where A_g is the geometric albedo in the Kepler bandpass. With this assumption in mind, and with the occultation depths of 21 ± 6 ppm and 22 ± 7 ppm for Kepler-5b and Kepler-6b, respectively (see Table 3), we find that both planets have low geometric albedos in the Kepler bandpass. Thus, we obtain A_g = 0.12 ± 0.04 and A_g = 0.11 ± 0.04 for Kepler-5b and Kepler-6b, respectively. Interestingly, we can infer the Bond albedo A_B from A_g. The Bond albedo is the fraction of the bolometric incident radiation that is scattered back out into space at all phase angle. Since we have no phase curve information for both objects, we assume a Lambertian criteria where A_B ⩽ 1.5 × A_g, which leads to A_B ⩽ 0.17 at the 1σ level for both object.

The equilibrium temperature is derived from

$$\begin{equation} T_{{\rm eq}}=T_\star (R_\star /a)^{1/2}[f(1-A_{\rm B})]^{1/4} \end{equation} \tag{ 1 }$$

which depends on the Bond albedo A_B and the re-distribution factor f which accounts for the efficiency of the transport of energy from the day to the nightside of the planet. f can vary between 1/4 for an extremely efficient redistribution (isothermal emission at every location of the planet) and higher values for an inefficient redistribution, implying big differences in the temperatures between the day and the nightsides of the planet. Nevertheless, the detection of the occultations in the Kepler bandpass may indicate that we measure the thermal emission of these planets at these wavelengths, as expected for the temperature regime of HJs. Therefore, the Bond albedo could be indeed well below 0.17. Notably, as Kepler spacecraft continues to monitor the photometry of Kepler-5b and Kepler-6b, the observations will provide better constraints on occultations and may reveal the phase curves which would allow to fully constrain the Bond albedo of these planets.

We estimate the thermal component of the planet's emission in the two Warm-Spitzer bandpasses from the occultation depths measured at 3.6 and 4.5 μm (see Table 2). To do so, we assume that the planetary emission is well reproduced by a blackbody spectrum and translate the measured depth of the secondary eclipse into brightness temperatures. We use the PHOENIX atmospheric code (Hauschildt et al. 1999) to produce theoretical stellar models for the star, with stellar temperature of T_eff = 6297 ± 60 for Kepler-5b (Koch et al. 2010b) and T_eff = 5647 ± 44 for Kepler-6b (Dunham et al. 2010). Taking the Warm-Spitzer spectral response function into account, the ratio of areas of the star and the planet and the stellar spectra, we derive the brightness temperatures that best fit the observed eclipse depths measured in the two IRAC bandpasses. The brightness temperature calculated this way resulted in T_Spitzer = 1930 ± 100 K and T_Spitzer = 1660 ± 120 K for Kepler-5b and Kepler-6b, respectively. If we further assume that the planet is in thermal equilibrium and has a zero Bond albedo, these temperatures favor low values of the re-distribution factor f, i.e., the planet dayside efficiently re-radiates the incoming stellar energy flux.

The brightness temperatures we derive using the Warm-Spitzer bandpasses are larger than the equilibrium temperatures for both targets (see Table 1). This may indicate that our assumption that these planets behave as blackbodies is most probably incorrect, similar to the giant planets of our solar system. Furthermore, a more robust determination of A_B requires detailed model computation because of the wide range of parameter space (Seager et al. 2005). We discuss below our results from two different types of atmospheric models.

In the first approach, we use the atmospheric retrieval method of Madhusudhan & Seager (2009) to derive the temperature structure and composition of each system, given the data. The model involves one-dimensional line-by-line radiative transfer, with parametric temperature structure and composition, and includes the major molecular and continuum opacity sources, along with constraints of LTE, hydrostatic equilibrium, and global energy balance. This modeling approach allows one to compute large ensembles of models (∼10⁶) and explore the parameter space of molecular compositions and temperature structure in search of the best-fitting models. In the present work, the dominant sources of opacity included are H₂O, CO, CH₄, CO₂, NH₃, and H₂–H₂ collision induced absorption in the infrared, and TiO, VO, Na, K, and Rayleigh scattering in the visible. In the present context, the number of model parameters is N = 10–15 (Madhusudhan & Seager 2009). Thus, the limited number of available observations (N_obs = 3) imply a substantial degeneracy in solutions. Consequently, our goal with the present data is to find a nominal set of solutions, as opposed to finding unique fits.

Models fitting the observations of Kepler-5b and Kepler-6b are shown in Figures 6 and 7. We find that observations of Kepler-5b can be explained to good precision by a wide range in composition, including models with equilibrium chemistry assuming solar abundances. However, current observations tentatively rule out a thermal inversion in this system. At the temperatures of Kepler-5b, the atmosphere is expected to be abundant in CO (Burrows & Sharp 1999), which has a strong feature in the 4.5 μm bandpass. A thermal inversion would therefore naturally predict a high flux in this channel over the 3.6 μm channel, as opposed to the observed fluxes which are similar between the two channels. Instead, the CO absorption feature caused by the lack of a thermal inversion explains the observations very well. Additionally, the Warm-Spitzer and Kepler data together require that the incident energy is mostly re-radiated on the dayside, with a low albedo and a small fraction (f) of dayside energy redistributed to the night side. Assuming zero albedo, the model shown has f = 0.12.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The observations of Kepler-6b allow for the possibility of a thermal inversion in its atmosphere. As shown in Figure 7, the data can be explained to within the 1σ errors by models with and without thermal inversions. As explained above, the higher flux in the 4.5 μm channel over the 3.6 μm channel can be interpreted as a sign of an inversion. However, given the degrees of freedom allowed by the molecular abundances, the limited data can be fit just as well without inversions (see, e.g., Madhusudhan & Seager 2011). A wide range of abundances can fit the data, including those close (within a factor of 10) to equilibrium chemistry assuming solar abundances. However, as with Kepler-5b, the Warm-Spitzer and Kepler data together place stringent constraints on the energy budget on the dayside. Our findings require low albedos and low redistribution to the nightside of this planet.

The second model we consider to compare our data to is the HJ atmospheric model of Fortney et al. (2008). Our motivation here is that, given our limited wavelength coverage, it is important to compare our best-fit planetary-to-star flux ratio values with physically motivated models. We aim at broadly distinguishing these atmospheres between different classes of models, given that there are no prior constraints on the basic composition and structure.

The atmospheric spectrum calculations are performed for one-dimensional atmospheric pressure–temperature (P–T) profiles and use the equilibrium chemical abundances, at solar metallicity, described in Lodders (2002) and Lodders & Fegley (2006). This is a self-consistent treatment of radiative transfer and chemical equilibrium of neutral species. The opacity database is described in Freedman et al. (2008). So far, this model has been used to generate P–T profiles for a variety of close-in planets (Fortney et al. 2005, 2006, 2008). The models are calculated for various dayside-to-nightside energy redistribution parameters (f) and allow for the presence of TiO at high altitude, playing the role of absorber, which likely leads to an inversion of the T–P profile.

We compare the data to model predictions and select models with the best reduced χ² . We note that this is not a fit involving adjustable parameters. All models for both targets show that the dayside is re-radiating the stellar flux efficiently, favoring low values for f. Since our observations are weakly constraining, model comparisons for Kepler-6b (Figure 8) suggest that both inverted and non-inverted atmospheric temperature profiles can reproduce the data. However, this is not the case for Kepler-5b where the models show no inversions. Nevertheless, the measurements in the visible are obtained by averaging approximately one year of Kepler data. During this time the planet emission and longitudinal temperature structure may change as suggested by several relevant theoretical studies (e.g., Cho et al. 2003; Showman et al. 2008; Rauscher & Menou 2010).

Zoom In Zoom Out Reset image size

For HJs such as Kepler-5 and Kepler-6, both scattered stellar light and planetary thermal emission could contribute to the planet emergent flux in the Kepler bandpass (López-Morales & Seager 2007). The ratio of scattered starlight to thermal emission depends on the atmospheric composition. These objects could be too hot for condensates. Without a reflective condensate layer, photons in the Kepler bandpass are absorbed before being scattered. Therefore, the albedo of Kepler-5 and Kepler-6 are likely to be low, as observed here.

Knutson et al. (2010) shows that there could be a correlation between the host star activity level and the thermal inversion of the planetary atmosphere. Effectively, the strong XEUV irradiation from the active stellar host of HJs could deplete the atmosphere of chemical species responsible for producing inversions. The Ca ii H & K line strengths are a good indicator of the stellar activity. We obtain Keck High Resolution Echelle Spectrometer spectra for Kepler-5 and Kepler-6, and estimate the line strengths of S_HK = 0.14 and log R'_HK = −4.97 for Kepler-5 and S_HK = 0.16 and log R'_HK = −4.98 for Kepler-6. Both stars are moderately quiet sub-giants, with log R'_HK falling in between cases for thermal inversion or no inversion (see discussion in Knutson et al. 2010). We also evaluate the empirical index defined by Knutson et al. (2010), which provides a way to distinguish between the different HJ atmospheres. Knutson et al. (2010) proposes an index value that could be correlated with the presence on a thermal inversion. Using the same definition, we find that the index = −0.023 ± 0.023 and index = 0.040 ± 0.033 for Kepler-5 and Kepler-6, respectively. These values suggest that the atmosphere of Kepler-6b is consistent with a weak thermal inversion and the one of Kepler-5b with a non-inverted profile, which are both in agreement with the results of our present study.

As the Kepler spacecraft continues to monitor the photometry of Kepler-5b and Kepler-6b, the future quarters of observations will be used to improve the signal-to-noise. This will provide better constraints on the occultations and may reveal the phase curves which would allow us to derive a more accurate Bond albedo and thus provide better constraints on the energy budget of these planetary atmospheres.