STUDYING ATMOSPHERE-DOMINATED HOT JUPITER KEPLER PHASE CURVES: EVIDENCE THAT INHOMOGENEOUS ATMOSPHERIC REFLECTION IS COMMON

We use here Kepler long-cadence data from quarters 2 through 17.⁵ Quarters 0 and 1 are ignored since they show a larger scatter than the other quarters. We have applied several preprocessing steps to the data, following Shporer et al. (2011, 2014), to prepare it for the atmospheric model fitting. We first remove instrumental signals, or trends, by fitting the first four cotrending basis vectors (CBVs) to the data of each quarter using the Pyke Python package (Still & Barclay 2012). We checked that increasing to six CBVs gave consistent results well within $1\sigma$, consistent with Shporer et al. (2011, 2014) and Hu et al. (2015).

We then fit each continuous data segment with a fifth-degree polynomial while ignoring in-transit and in-occultation data, and then we divide all data within the segment by that polynomial. This step does not affect the sinusoidal modulation along the orbit since the duration of each continuous segment is at least an order of magnitude longer than the orbital period. Fitting was done while iteratively rejecting $5\sigma$ outliers until none are identified using the formalism described in Shporer et al. (2014, see their Section 2). Using a polynomial of degrees 4 and 6 did not change the results. Using smaller degrees showed an increased overall scatter in the light curves while showing the same sinusoidal amplitude, where the latter is the astrophysical signal we seek to characterize here. Using larger degrees decreased the sinusoidal amplitude resulting from the excessively aggressive detrending. The phase curves before and after the polynomial detrending are presented in Appendix A.

Next, we phase folded each light curve using the known ephemeris and binned it using 200 phase bins. We rejected $4\sigma$ outliers from each bin, which contained about 230 and 300 individual Kepler measurements for Kepler-12 and Kepler-41, respectively. For each bin, we calculated the median relative flux value and took the median absolute deviation (also known as MAD) times 1.4826 to be the uncertainty. Therefore the bins' uncertainties are independent of the error bars assigned to the Kepler individual measurements and are based on the true scatter in each phase bin. We compared the measured scatter within each phase bin to the expected scatter based on the error bars of the individual measurements and Poisson statistics. We found that on average the measured scatter is 11% and 49% larger than the expected scatter for Kepler-12 and Kepler-41, respectively. For Kepler-41, the increased measured scatter can be at least partially attributed to stellar activity (see Section 4.3). Although, for Kepler-12 there is no evidence of stellar activity (see Section 4.3), so the increased measured scatter is attributed to residual correlated noise, which in principle may also exist in the Kepler-41 data, but for the latter the stellar activity signal dominates.

Next, we checked whether combining the data from different Kepler quarters follows Poisson statistics or whether it adds some noise due to, e.g., instrumental differences between quarters or data detrending. For each object, we generated a phase-folded and binned light curve for each quarter. We then calculated for each phase bin the ratio between the average uncertainty, averaging across all quarters, and the uncertainty of the phase bin when using all Kepler data, and further dividing by the square root of the number of quarters. For both the Kepler-12 and Kepler-41 data sets, the average ratio across all phase bins was within 1%–2% from unity, showing that combining data from different quarters closely follows Poisson statistics.

Finally, we checked that our results do not depend on the phase bin width by also analyzing phased light curves with 100 and 400 bins and confirming the results are identical.

Before moving to the atmospheric modeling, we removed from the phase-folded light curves the orbital modulations induced by beaming and ellipsoidal, based on the known planet mass in each system (Demory et al. 2011; Santerne et al. 2011). This is done for completeness since, as mentioned above, for the systems analyzed here these modulations are small, at the noise level, and the phase-curve amplitude is dominated by atmospheric processes (see Figure 1), so any approximations in determining the shape of the beaming and ellipsoidal phase curves do not affect the overall phase-curve shape significantly.

Zoom In Zoom Out Reset image size

We analyze the Kepler phase curves using the interpretation framework described in Hu et al. (2015). The model efficiently computes reflection and thermal emission contributions to the phase curve, considering both hot spot shift due to equatorial winds and inhomogeneous atmospheric reflection due to patchy clouds. The model assumes the atmosphere has clear longitudes and cloudy longitudes, and the cloudy part is more reflective than the clear part. The model assumes the cloud distribution is controlled by the temperature distribution, which is in turn approximated by the analytical model proposed by Cowan & Agol (2011). The model also assumes the thermal emission comes from a thermal photosphere that may be hotter than the equilibrium temperature. Here, the equilibrium temperature, as a function of longitude, is defined to be the atmospheric temperature without any greenhouse effects and is calculated from the Bond albedo and the heat redistribution efficiency.

We fit the phase curve to this model characterized by five parameters: the Bond albedo, a heat redistribution efficiency, a greenhouse factor, the condensation temperature of clouds, and a reflectivity boosting factor by clouds. The Bond albedo describes the overall reflectivity of the planet and takes any value between 0 and 1. The heat redistribution efficiency, the ratio between the radiative timescale and the advective timescale, controls the longitudinal temperature distribution. The magnitude of this parameter describes how well heat is transported: the larger the parameter, the smaller the longitudinal temperature variation becomes. The sign of this parameter describes the direction of the equatorial winds: a positive sign indicates super-rotating eastward winds, and a negative sign indicates westward winds. The greenhouse factor is the ratio between the temperature of the thermal photosphere and the equilibrium temperature and takes any value not smaller than 1. The cloud condensation temperature, another free parameter, determines the longitudinal boundaries of the clouds. Finally, the reflectivity boosting factor by clouds is the proportional increase in reflectivity of an atmospheric patch when it becomes cloudy and takes any positive value. The rest of the system parameters are known and taken from the literature (Fortney et al. 2011; Santerne et al. 2011; Quintana et al. 2013).

We use the Markov Chain Monte Carlo (MCMC) method (Haario et al. 2006) to explore the parameter space. Following Gelman & Rubin (1992), we calculate two Markov chains for each of the two objects, each chain containing 10⁶ steps where the first half of each chain is considered the "burn-in" period and is removed once the chain is completed. We verify that the R values for all parameters are less than 1.01 to ensure convergence (Gelman & Rubin 1992), and that in both cases the results from the two chains are identical. The allowed ranges and the prior distributions of the model parameters are as follows: the Bond albedo uniformly ranges in [0, 1]; the heat redistribution efficiency uniformly ranges in [0, 100]; the greenhouse factor uniformly ranges in [1, 2]; the cloud condensation temperature uniformly ranges in [1000, 3000] K; and the reflection boosting factor uniformly ranges in [0, 100]. The redistribution efficiency parameter is only allowed to take a positive value because atmospheric circulation models predict eastward equatorial winds (e.g., Showman & Guillot 2002; Showman & Kaspi 2013). This choice effectively requires any post-occultation phase offset to be explained by the asymmetric reflection components induced by the patchy clouds, rather than the thermal emission component (Hu et al. 2015). Based on the converged Markov chain, we derive the occultation depth, the phase-curve amplitude, and the phase shift of the phase-curve maximum from occultation phase. The latter is defined to be positive for post-occultation maximum.

After carrying out the analysis once, we noticed that for Kepler-12b the reduced ${\chi }^{2}$ of the residuals was 1.12. Therefore, we repeated the analysis of that phase curve while increasing the bins' uncertainties to bring the reduced ${\chi }^{2}$ to unity. This has increased the uncertainty on the fitted parameters by up to a few percent. For Kepler-41b the reduced ${\chi }^{2}$ was 0.86 in the first analysis, so no further analysis was done.

To test our results, we applied the so-called "prayer bead" approach to both objects, where the residuals are cyclicly permuted and added back to the model, and the new data set is then refitted. This approach preserves any correlated noise features in the data and produces a distribution for each fitted parameter from which its uncertainty can be estimated (e.g., Bouchy et al. 2005; Southworth 2008; Winn et al. 2008). For both objects, the prayer bead analysis did not result in larger uncertainties than the original MCMC analysis, so we use the latter as our final result.

We measure a phase-curve amplitude of 17.5 ± 2.2 ppm (parts per million) and an eclipse depth of 14.3 ± 2.5 ppm, consistent with previous studies (Angerhausen et al. 2014; Esteves et al. 2015) but with smaller uncertainties. We detect a clear phase-curve asymmetry, with a phase-curve maximum shifted by $47\buildrel{\circ}\over{.} {6}_{-7.6}^{+6.9}$ from the occultation phase. This is consistent with the phase shift reported by Esteves et al. (2015) but $\approx 2\sigma$ smaller than that of Angerhausen et al. (2014; D. Angerhausen 2015, private communication). We find that the phase-curve amplitude is greater than the eclipse depth by ≈$1\sigma$, in agreement with the measurements of Esteves et al. (2015). The difference between the phase-curve amplitude and eclipse depth is enhanced by the relatively large phase offset.

To further test our analysis, we carried out a separate analysis where we fitted a simple parametric model including a sinusoidal amplitude and phase. This is similar to previous analyses (see a summary of previous results in Table 9 of Esteves et al. 2015). We derived a peak-to-peak amplitude of 19.0 ± 2.2 ppm and a phase of maximum shifted later than the occultation phase by 58.5 ± 6.8 deg. These values are ≲$1\sigma$ away from that of our original analysis and with similar uncertainties (see Table 1).

Looking at the Kepler-12b phase curve (Figure 1 left panel), there are a few consecutive bins around phase 0.2 that fall below the fitted model. While it could be an astrophysical signal that is not included in our model, it could also be a statistical, correlated noise feature. Although, correlated noise is expected to average out in the phase-folded light curve. Either way, as described in Section 2.2, we have inflated the bins' uncertainties to bring the residuals reduced ${\chi }^{2}$ to unity, and the prayer bead analysis did not result in larger error bars for the fitted parameters compared to the MCMC analysis.

Although an asymmetry in the Kepler-41b phase curve is not identified by Angerhausen et al. (2014), it is identified by Esteves et al. (2015) and is visually apparent in Santerne et al. (2011, their Figure 3) and Quintana et al. (2013, their Figure 2). Here, we confirm the post-occultation phase maximum that we measure to be shifted by 178 ± 63. We derive an eclipse depth of ${49.1}_{-5.2}^{+6.4}$ ppm and a phase amplitude of 48.7 ± 5.2 ppm, within 1–2σ from previous results but with typically smaller error bars (Santerne et al. 2011; Quintana et al. 2013; Angerhausen et al. 2014; Esteves et al. 2015). Esteves et al. (2015) report a phase amplitude greater than the eclipse depth (although at low statistical significance of only $1.2\sigma$), while here we find a phase-curve amplitude that is highly consistent (within $0.1\sigma$) with the eclipse depth, indicating either poor heat redistribution or large planetary albedo.

We fitted Kepler-41b with the same simple parametric model as we did for Kepler-12b, consisting of a sinusoidal amplitude and phase. We derived a peak-to-peak sinusoidal amplitude of 39.2 ± 4.2 ppm and a phase of maximum shifted later than the occultation phase by 252 ± 59. These are within $1.5\sigma$ from the results of our original analysis and with similar uncertainties (see Table 1).

Our analysis of Kepler-12b and Kepler-41b Kepler phased light curves confirms and refines previous results and includes a more detailed atmospheric analysis than done previously.

We identify an asymmetry in the Kepler phased light curves for the two planetary systems studied here, where the maximum is post-occultation. Our interpretation is that this asymmetry is due to the brightest region on the planetary surface being shifted westward from the substellar point that in turn originates from an asymmetric surface distribution of clouds in the planetary atmosphere. The post-occultation maximum requires that the reflective clouds exist on the west side of the substellar point but not on the east side. Such asymmetric cloud distribution may be a result of eastward heat transport driven by zonal jets in the atmosphere. In addition, the clear part of the atmosphere must be quite dark, probably due to absorption by Alkali metal vapors deep in the atmosphere, whereas the cloudy part appears to be bright due to Mie scattering (Sudarsky et al. 2000, 2003).

For Kepler-41b, in order to produce a phase variation consistent with Kepler observations, the east side of the substellar meridian should be devoid of reflective clouds and have a reflectivity less than 0.2 ($1\sigma$ limit; see Table 1). Figure 1 shows that the thermal emission component has the opposite phase offset than the asymmetric reflection component. On the west side, where reflective clouds are present, the reflectivity should be greater than 0.3. The overall Bond albedo of the planet is constrained to be 0.25 ± 0.06. The cloud condensation temperature parameter is found to be ${1580}_{-120}^{+150}$ K. To compare with material properties, one should also consider that clouds can form deep in the atmosphere and hence at a higher temperature (Hu & Seager 2014). After multiplying by the derived greenhouse factor, the $1\sigma$ cloud condensation temperature range is 1530–1910 K. Silicates (MgSiO₃ and Mg₂SiO₄) are consistent with this temperature range and are highly reflective, making them candidate condensates.

Comparing Kepler-12b with Kepler-41b, we find that Kepler-12b is a more extreme case. Its phase curve is completely dominated by asymmetric surface distribution of reflective clouds, whereas the phase curve of Kepler-41b is dominated by symmetric reflection (see Figure 1). Because the cloud distribution is controlled by the cloud condensation temperature, the post-occultation phase offset yields a tight constraint on the latter (${1506}_{-24}^{+13}$ K). Multiplying by the derived greenhouse factor, the $1\sigma$ cloud condensation temperature range is 1480–1670 K. To produce the phase offset, the clear part of Kepler-12b atmosphere needs to have an albedo very close to zero, lower than that of Kepler-41b. This can be naturally explained as a higher metallicity and therefore greater Na and K absorption in the atmosphere of Kepler-12b (Sudarsky et al. 2000, 2003). However, the host star in the Kepler-12 system is relatively metal poor compared to Kepler-41 (Fortney et al. 2011; Santerne et al. 2011); therefore, it is unclear why Kepler-12b would have a more metal-rich atmosphere than Kepler-41b.

Overall, our analysis shows that it is reasonable to assume that an asymmetric cloud distribution is a common phenomena for hot Jupiters, but the details for the clear and cloudy parts of their atmospheres can still manifest great diversities.

For both planetary atmospheres, the current data constrain the Bond albedo and the cloud condensation temperature and put an upper limit on the greenhouse factor. Among these parameters, we see that the Bond albedo and the cloud condensation temperature are somewhat correlated, as expected (see Figures 4 and 5). This is because to produce a patchy cloud solution, the cloud condensation temperature needs to be between the maximum and the minimum temperature of the dayside, both affected by the Bond albedo. The same correlation is observed for Kepler-7b (Hu et al. 2015). We note that the efficiency of heat redistribution cannot be sufficiently constrained by the current data. Although the best-fit models suggest efficient heat redistribution, many more models would be permitted at an expense of minimally increasing ${\chi }^{2}$. For example, we find that much less efficient heat redistribution and a much greater reflectivity contrast between the cloudy part and the clear part could render an almost-as-good fit for Kepler-41b. As discussed in Hu et al. (2015), this insensitivity to the exact value of the heat redistribution efficiency is tied to the uncertainty in the exact condensation temperature of the cloud material.

The atmospheric analysis done here for Kepler-12b and Kepler-41b and the analysis done for Kepler-7b by others (Demory et al. 2013; Hu et al. 2015; Webber et al. 2015) point to a similar scenario: that the optical Kepler phase curve is dominated by reflected light (as opposed to thermal emission) and that the shift between the phase-curve maximum and planetary occultation phase is a result of a non-uniform longitudinal cloud coverage. The irradiation received by these three planets is similar, and hence the cloud condensation temperatures derived from our phase-curve analysis are also similar. To form the patchy cloud prominent longitudinal distribution, the condensation/evaporation timescale of the cloud particles needs to be shorter than their advection timescale (Hu et al. 2015). This condition is uniquely met by hot Jupiters that have fast zonal winds in their atmospheres (e.g., Knutson et al. 2007).

The three planets identified here as having a phase curve dominated by atmospheric processes, hence they allow a direct view of their atmospheres, all show a phase shift in the form of a post-occultation maximum. The identification of phase shifts in these atmosphere-dominated phase curves suggests they are common and exist also in phase curves of other hot Jupiters where they are more difficult to identify uniquely as those systems show additional phase modulations induced by beaming and tidal ellipsoidal distortion. This has several potential implications that we discuss below.

One obvious implication of the phase shift is on atmospheric studies using only the occultation and not the entire phase curve. Such studies commonly assume that the occultation depth is a flux measurement of the brightest region on the planetary surface, which is not the case for shifted phase curves. In those cases, the secondary eclipse depth underestimates the planetary brightest-region flux, which in turn can lead to underestimating the planetary reflectivity and/or thermal emission.

Another implication of the phase shift involves the BEER model (Faigler & Mazeh 2011), developed to search for non-transiting systems through the detection of photometric orbital phase modulations. It assumes that all processes, including gravitational (beaming and ellipsoidal) and atmospheric (reflected light and thermal emission), are all well aligned with the phase of transit and occultation and that the orbit is circular. Meaning, the extrema of each effect occurs exactly at either transit, occultation, or quadrature phases. When the reflection component is shifted, with either a post- or pre-occultation maximum, then the BEER model will give a biased estimate of the beaming amplitude. This is so since the beaming and reflection components are the sine and cosine components of the same frequency modulation (see Faigler & Mazeh 2011 for more details). This in turn will make the companion mass estimate based on the beaming amplitude, m_2,beam, disagree with the estimate based on the ellipsoidal amplitude, m_2,ellip. A post-occultation maximum, as identified here, will decrease the beaming amplitude measured by BEER, leading to m_{2, beam} < m_{2, ellip}, while a pre-occultation maximum will lead to m_{2, beam} > m_{2, ellip}.

As noted earlier, the mass discrepancy between the beaming-based mass estimate and the ellipsoidal-based mass estimate has already been measured for several systems, including transiting planets and stellar binaries with a transiting white dwarf companion (see Section 1). Another source of disagreement between the beaming-based and ellipsoidal-based mass estimates are the approximations in the interpretation of the ellipsoidal light-curve shape (Morris 1985; Morris & Naftilan 1993; Pfahl et al. 2008), especially in cases where the primary star rotation is not synchronized with the orbit or its spin axis is not aligned with the orbital angular momentum axis. Other factors that might be further contributing to this discrepancy are poor understating of the host star (e.g., Shporer et al. 2014) and non-zero orbital eccentricity. Therefore, deciphering the origin of the discrepancy requires a careful case-by-case study. Faigler & Mazeh (2015) recently suggested that the discrepancy for Kepler-13Ab, TrES-2b, HAT-P-7b, and Kepler-76b identified in Kepler data can be attributed to a pre-occultation phase-curve maximum of the atmospheric component. They use the difference between the two mass estimates to derive the angular shift of the brightest region on the planetary surface from the substellar point and attribute this shift to asymmetric thermal emission of the planetary atmosphere.

Yet another possible implication of the phase shift and the planetary asymmetric surface brightness distribution is an asymmetry, or a distortion, of the occultation light-curve ingress and egress. When identified at sufficiently high signal-to-noise ratio, it can be used to map the planet's surface brightness, a method known as eclipse mapping (e.g., Williams et al. 2006; de Wit et al. 2012). When the distortion is not directly identified, it can still lead to a shift in the measured occultation phase when modeled using a symmetric occultation light-curve model (e.g., Shporer et al. 2014). However, the Kepler data of the three transiting planets studied here are not sensitive to this minute effect.

Finally, we note that a longitudinally non-uniform atmosphere will affect the shape of the polarization reflected light phase curve (e.g., Seager et al. 2000; Berdyugina et al. 2011; Madhusudhan & Burrows 2012).

The three exoplanets studied here (Kepler-7b, Kepler-12b, and Kepler-41b) all have a mass of roughly half a Jupiter mass and orbital semimajor axes within 0.1 AU. Such a planet mass gives a gravitational signal (due to the beaming and ellipsoidal effects) at the ∼1 ppm level, which is not expected to be detectable in Kepler data. On the other hand, their Kepler light curves show a clear atmospheric photometric orbital modulation, which, if detectable, can reveal the existence of a faint low-mass companion even when the system does not transit, as discussed in greater detail below. Therefore, such planets represent a class of non-transiting planets that although are typically beyond the reach of the BEER approach (Faigler & Mazeh 2011), as the latter is designed to simultaneously detect all orbital modulations, can potentially still be identified through their atmospheric photometric orbital modulation resulting from their relatively large geometric albedo. The latter is similar to the approach discussed in detail by Jenkins & Doyle (2003).

Since the overall phased light-curve shape is similar to a pure sine, a simple Fourier-like period analysis could be used to identify it, at least in principle. This will result in a measurement of the orbital period, and the photometric modulation amplitude will give a mutual constraint on the companion's radius, geometric albedo, and brightness temperature at the observed wavelength. The companion's (minimum) mass and the true nature of the system can then be obtained via radial velocity follow-up measurements. Interestingly, Charpinet et al. (2011) and Silvotti et al. (2014) used a similar approach to claim the discovery of non-transiting small planets orbiting the hot B subdwarf stars Kepler-70⁶ and Kepler-429,⁷ respectively, at short orbital periods.

To see whether Kepler-7b, Kepler-12b, and Kepler-41b could have been identified as candidate hot Jupiter systems, even if they were not transiting, we show in Figure 2 a Lomb–Scargle (L–S; Lomb 1976; Scargle 1982) periodogram of their Kepler light curves after removing the transit and occultation data. We used the fast L–S algorithm described by Press & Rybicki (1989) as implemented in the Matlab/plomb function. In all three panels, the dashed horizontal black line marks the L–S power level corresponding to a false-alarm probability of ${10}^{-6}$, and the frequency component at the orbital frequency is marked by a gray arrow.

Zoom In Zoom Out Reset image size

Kepler-41 (Figure 2, bottom panel) is known as an active star (Quintana et al. 2013), and its periodogram is dominated by the rotational modulation at frequencies smaller than the orbital frequency. Yet, there is a clear signal at the orbital frequency. For Kepler-7 and Kepler-12 (Figure 2, top and middle panels), the strongest frequency components are at the orbital frequency. Meaning, they would have been identified as periodic variables at the orbital period even with no transits. In fact, including the data during conjunctions, which were removed here, would make the orbital frequency component even stronger, as we have confirmed through simulations.

However, there are other processes that may result in photometric sinusoidal variability, so the simplistic approach described above may suffer from a high false-positive rate (Jenkins & Doyle 2003). One such scenario involves an active star, where the combination of the stellar rotation with the non-uniform spot surface distribution results in sine-like photometric modulations (e.g., Hartman et al. 2011; McQuillan et al. 2014). Although, such modulations vary with time in amplitude and phase, as opposed to a steady orbital signal, so a more sophisticated approach, e.g., using wavelets (Bravo et al. 2014), can differ between activity-induced modulations and orbital modulations. Moreover, as shown by McQuillan et al. (2013), hot Jupiters orbiting at short periods are typically hosted by stars rotating at significantly longer rotation periods, as is the case for Kepler-41b. So the orbital and rotational frequency components are expected to be well separated in the periodogram.

Stellar asteroseismic pulsations might also be a false-positive scenario, although, typical pulsations of main-sequence stars are at a much shorter timescale than hot Jupiter orbital periods (e.g., Chaplin & Miglio 2013).

Yet another false-positive scenario can originate from blending of nearby stars on Kepler pixels. In that scenario, the light from a faint variable star, for example, an ellipsoidal variable (contact stellar binary), can be blended with that of a brighter non-variable star to give an observed low-amplitude sinusoidal modulation. This scenario, which presents a problem for any variability study using Kepler data, can be examined using light curves of the individual Kepler pixels (Bryson et al. 2013) combined with high angular resolution imaging of the target that maps the stars within the Kepler pixel mask (e.g., Horch et al. 2014; Law et al. 2014).

Although the detection of hot Jupiters through their atmospheric modulations alone is expected to be more prone to false-positive scenarios than BEER, it can allow access to (non-transiting) exoplanets in a lower-mass range, down to below a Jupiter mass, where exoplanets are far more common than at a few Jupiter masses or above where they become detectable by the BEER approach (Shporer et al. 2011; Faigler et al. 2013).