© ESO, 2016

1. Introduction

Many observations reported in recent years are compatible with the presence of clouds in exoplanetary atmospheres.

Transmission spectra measure the transit depth during primary transit as a function of wavelength. If the atmosphere is opaque at a given wavelength, the planet will appear larger at that wavelength. However, visible, near-IR and IR transmission spectra of planets such as GJ1214b (e.g., Bean et al. 2010; Kreidberg et al. 2014), GJ436b (e.g., Knutson et al. 2014), CoRoT-1b (e.g., Schlawin et al. 2014), WASP-31b (e.g., Sing et al. 2015), or HAT-P-32b (e.g., Gibson et al. 2013), show very little variation with wavelength. They can be approximated by flat lines to within measurement uncertainties. This has been interpreted as being due to optically thick cloud layers high up in the atmosphere that mask the expected molecular and atomic absorption.

During secondary eclipse, the apparent dayside brightness temperature can be measured at IR wavelengths. Again, different wavelengths probe different atmospheric pressures due to wavelength-dependent opacity, and so brightness temperatures are expected to vary. As was observed for transmission spectra, IR emission spectra of many exoplanets are consistent with being featureless (e.g., Hansen et al. 2014). A single blackbody brightness temperature can be assigned to the photosphere, independent of wavelength. This is consistent with uniform cloud coverage that blocks the access to deeper atmospheric layers.

UV-visible secondary eclipse spectra have been used to determine the wavelength-dependent albedo of hot Jupiters (e.g., Evans et al. 2013). It is found that the slope in the albedo spectrum is consistent with clouds obscuring deeper atmosphere levels where molecular and atomic absorption could influence the observations.

The CoRoT and Kepler space missions have provided a wealth of high-precision, broadband optical photometry of thousands of exoplanets and candidates. About 20 of these exoplanets show detectable phase curves (e.g., Snellen et al. 2009; Mazeh & Faigler 2010; Barclay et al. 2012; Quintana et al. 2013; Esteves et al. 2013, 2015), i.e., variations with orbital phase as the planetary dayside (which reflects/re-emits starlight back to the observer) rotates in and out of view. A few Kepler phase curves show asymmetries with respect to secondary eclipse. Post-eclipse maxima of the phase curves have been found for at least four exoplanets (e.g., Demory et al. 2013; Esteves et al. 2015; Webber et al. 2015). These are conceptually interpreted as being a consequence of inhomogeneous cloud cover on the planet dayside. Simply put, clouds form on the nightside and evaporate as they are transported from the cooler “morning” side across the substellar meridian to the hotter “evening” side¹.

In recent years, detailed theoretical cloud modeling has been performed (e.g., Sudarsky et al. 2000; Parmentier et al. 2013; Wakeford & Sing 2015; Webber et al. 2015; Lee et al. 2015). These studies used temperature structures calculated by 1D and 3D atmospheric models to predict cloud condensate composition (e.g., silicate and iron clouds have been proposed) and cloud altitude (usually around 1–100 mbar of pressure). Together with 3D circulation patterns, sophisticated cloud formation schemes simulate sedimentation and cloud particle growth to calculate cloud particle size distributions, cloud particle densities and cloud location. These studies generally confirm the heterogeneous cloud cover inferred from optical phase curves. However, given current data quality, degeneracies exist between cloud parameters that are not easily remedied with optical phase curves alone. For instance, even homogeneous cloud cover (or clouds clustered towards the substellar point) could produce asymmetric phase curves if the scattering asymmetry factor of the cloud particles is adapted accordingly.

Therefore, apart from intensive numerical modeling efforts, additional diagnostics for cloud properties are needed to help resolve the degeneracies and better constrain cloud characteristics. In this work, we propose a simple concept to probe cloud locations, i.e., asymmetric primary transits². The concept is based on the premise that, in the case of heterogeneous cloud cover, the western limb of the planet (the leading limb in primary transit) is cloudy, whereas the eastern limb (trailing limb) is cloud-free, as proposed based on phase curve observations. Thus, the planet is composed of two hemispheres with different apparent radii, which produces different transit signatures during ingress and egress.

Transit observations, especially from space-based surveys, are particular useful for this method for several reasons. Mainly, the quasi-continuous coverage allows for the analysis of hundreds of near-consecutive transits, hence an excellent time resolution of the order of seconds or less in the phase-folded light curves. Furthermore, compared to spectroscopic observations, broadband photometry produces much higher signal-to-noise ratios and, thus enables stronger constraints on cloud parameters.

We apply the new technique to two planets in the Kepler field with well-characterized asymmetric phase curves, namely Kepler-7b (e.g., Latham et al. 2010; Demory et al. 2011, 2013; Esteves et al. 2015) and HAT-P-7b (Pál et al. 2008; Borucki et al. 2009; Welsh et al. 2010; Esteves et al. 2015; von Paris et al. 2016). Furthermore, we analyze transit light curves of HD 209458b taken from Brown et al. (2001; their Fig. 3) obtained from four transits observed with the Hubble Space Telescope. HD 209458 is one of the brightest stars around which a transiting planet has been found, and the transmission spectrum of HD 209458b shows some evidence of possible cloud and/or haze cover (e.g., Deming et al. 2013). Therefore, we include it here for illustration purposes.

Although the results are not conclusive with respect to the presence of clouds, we will establish the applicability of the technique for future photometric surveys.

The paper is organized as follows: Sect. 2 presents the data reduction process adapted for this work, and Sect. 3 presents the numerical model used in the transit modeling. Section 4 introduces the cloud diagnostics for the primary transit light curve. The inverse model used to fit the forward models to the data is described in Sect. 5. Results are presented in Sect. 6 and discussed in Sect. 7. We conclude with Sect. 8.

2. Data reduction

For target planets observed by Kepler (i.e., Kepler-7b and HAT-P-7b), we use the short-cadence (SC) single-aperture photometry (SAP), which is publicly available from the MAST archive³. An example is shown in Fig. 1.

First, data points that are flagged by the Kepler pipeline as NaNs and obviously flawed data, are removed manually to produce a raw light curve R_qk for each quarter k. To obtain phase-folded transit light curves L_qk, we then proceed in three steps.

• 1.

  Outlier removal: We calculate a running median M_2.4 for each point j of the light curve, using a 2.4-h window, ${\mathit{M}}_{\mathrm{2.4}}\mathrm{\left(}\mathit{j}\mathrm{\right)}\mathrm{=}\mathrm{median}\left(\mathit{R}\mathrm{\left(}\mathit{i}\mathrm{\right)}\mathit{,}\mathrm{\left|}\mathit{t}\mathrm{\right(}\mathit{i}\mathrm{)}\mathrm{-}\mathit{t}\mathrm{(}\mathit{j}\mathrm{\left)}\mathrm{\right|}⩽\mathrm{1.2}\mathrm{h}\right)\mathit{,}$(1)where t denotes the time stamps of the data point. From there, we calculate the median absolute deviation (MAD) as follows: $\mathrm{MAD}\mathrm{=}\mathrm{median}\left(\mathrm{|}\mathit{R}\mathrm{-}{\mathit{M}}_{\mathrm{2.4}}\mathrm{|}\right)\mathit{.}$(2)A point j is rejected if $\mathrm{\left|}\mathit{R}\mathrm{\right(}\mathit{j}\mathrm{)}\mathrm{-}{\mathit{M}}_{\mathrm{2.4}}\mathrm{(}\mathit{j}\mathrm{\left)}\mathrm{\right|}⩾\mathrm{4}\mathrm{\times }\mathrm{MAD}\mathit{.}$(3)Figure 2 illustrates this procedure, for the 9th transit of Kepler-7b in Q3 (green lines).

  Fig. 2 Ninth transit in Q3: running 2-day median (red line) and 2.4-h median (green plain line). Green dashed lines encompass the 4xMAD criterion to identify outliers (see text).

• 2.

  Normalization: After outlier removal, we calculate a running median M₄₈ for each point j, using a roughly 2-day window (red line in Fig. 2): ${\mathit{M}}_{\mathrm{48}}\mathrm{\left(}\mathit{j}\mathrm{\right)}\mathrm{=}\mathrm{median}\left(\mathit{R}\mathrm{\left(}\mathit{i}\mathrm{\right)}\mathit{,}\mathrm{\left|}\mathit{t}\mathrm{\right(}\mathit{i}\mathrm{)}\mathrm{-}\mathit{t}\mathrm{(}\mathit{j}\mathrm{\left)}\mathrm{\right|}⩽\mathrm{24}\mathrm{h}\right)\mathit{.}$(4)The normalized light curve N is obtained by dividing the raw light curve by the median-filtered light curve: $\mathit{N}\mathrm{=}\frac{\mathit{R}}{{\mathit{M}}_{\mathrm{48}}}\mathrm{·}$(5)In this way, slow variations (planetary phase curve, stellar activity, etc.) are mostly removed from the final light curve. Figure 3 shows the normalized light curve for Kepler-7 in quarter 3.

  Fig. 3 Normalized SC-SAP for Kepler-7 in Q3.

• 3.

  Phase folding: the normalized light curve N is then phase-folded to obtain the final light curve L as a function of orbital phase Φ. For this, Φ is defined such that Φ = 0 at primary transit, $\mathrm{\Phi }\mathrm{\left(}\mathit{j}\mathrm{\right)}\mathrm{=}\frac{\mathit{t}\mathrm{\left(}\mathit{j}\mathrm{\right)}\mathrm{-}{\mathit{T}}_{\mathrm{0}}}{\mathit{P}}\mathrm{-}\mathrm{\lfloor }\frac{\mathit{t}\mathrm{\left(}\mathit{j}\mathrm{\right)}\mathrm{-}{\mathit{T}}_{\mathrm{0}}}{\mathit{P}}\mathrm{\rfloor }\mathit{,}$(6)where T₀ is the mid-transit time of the first transit, P the orbital period and ⌊ x ⌋ represents the floor function, i.e., the greatest integer less than or equal to x. Figure 4 shows the final light curve of Kepler-7b, using all quarters with SC data.

The values for T₀ and P are taken from the latest data release (DR24)⁴ and summarized in Table 1.

Figure 5 shows the Q0-Q17 light curve of HAT-P-7b. Since HAT-P-7 is about 2.5 mag brighter than Kepler-7 in the Kepler bandpass, the resulting light curve is much cleaner.

Figure 6 shows the transit light curve of HD 209458b, used in this work. Data is taken from Brown et al. (2001; their Fig. 3). Observations were made with the Hubble Space Telescope during four transits of HD 209458b and then phase-folded to obtain a composite light curve. As is clearly seen, the photometric quality of the light curve is exceptional, since HD 209458 is a V = 8^m star and was observed with the HST (aperture 2.4 m). However, an important point is the relatively low time resolution compared to the Kepler light curves.

3. Forward model

In this model, the planet is assumed to be composed of two hemispheres with radius R_trail and R_lead for the trailing and leading hemisphere, respectively. These two radii are parameterized as follows: $\begin{array}{ccc}{\mathit{R}}_{\mathrm{trail}}& \mathrm{=}& {\mathit{R}}_{\mathrm{p}}\\ {\mathit{R}}_{\mathrm{lead}}& \mathrm{=}& {\mathit{R}}_{\mathrm{p}}\mathrm{+}{\mathit{h}}_{\mathrm{cloud}}\mathit{,}\end{array}$where R_p is the planetary radius and h_cloud is the cloud altitude. We note that h_cloud< 0 is permitted, which then implies a larger trailing hemisphere. The concept of the model is illustrated in Fig. 7.

To calculate the stellar total luminosity and estimate the area obscured by the planet, we use a N_∗ × N_∗ grid to resolve the stellar disk, with its center corresponding to the stellar center (see Fig. 7). N_∗ is chosen such that a stellar pixel corresponds to roughly half a scale height of the planet in question (N_∗ = 4000 for Kepler-7b, N_∗ = 10 000 for HAT-P-7b, for example).

As is general practice in exoplanet transit modeling, the asymmetric model adopts a quadratic limb-darkening law with u₁ and u₂ as the linear and the quadratic limb-darkening coefficients for the intensity distribution I across the stellar disk: $\mathit{I}\mathrm{\left(}\mathit{\mu }\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{-}{\mathit{u}}_{\mathrm{1}}\mathrm{·}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\mu }\mathrm{)}\mathrm{-}{\mathit{u}}_{\mathrm{2}}\mathrm{·}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\mu }{\mathrm{)}}^{\mathrm{2}}\mathit{,}$(9)with μ = cos(φ) and φ is the angle between the local surface normal and the observer direction.

Planetary orbits are assumed to be circular, which is a good approximation for the targets considered (see, e.g., Winn et al. 2009; Demory et al. 2011, 2013; Esteves et al. 2015).

Further parameters needed to describe the passage of the planet in front of the star are the impact parameter b and a_S, the projected semi-major axis a_S in units of the stellar radius: ${\mathit{a}}_{\mathrm{S}}\mathrm{=}\frac{\mathit{a}}{{\mathit{R}}_{\mathrm{\ast }}}\mathit{,}$(10)where a is the semi-major axis and R_∗ is the stellar radius (fixed in this work due to strong asteroseismology constraints, e.g. Demory et al. 2011; Van Eylen et al. 2012).

The impact parameter is defined as the projected distance (in units of stellar radii) of the planets’ center to the stellar center at mid-transit. It is related to the orbital inclination i via: $\mathit{b}\mathrm{=}\cos \mathrm{(}{\mathit{i}}^{\mathrm{)}}\mathrm{·}{\mathit{a}}_{\mathrm{S}}\mathit{.}$(11)In addition, we introduce a seventh parameter which is the phase offset φ₀ of the passage of the planetary meridian over the stellar meridian (see Fig. 7) compared to the minimum of the primary transit. In general, we assume φ₀ = 0, however due to, e.g., eccentric orbits, φ₀ can be different from zero.

The passage of the planet in front of the star is discretized such that the planet advances less than one stellar pixel per time step.

We have compared the new model with the widely used, analytical Mandel & Agol (2002) model. Usually, the deviations between both models are of the order of 1 ppm or less, hence much smaller than the signals we look for (see next section, Fig. 12). We deem such a deviation acceptable. With increasing numerical resolution in our model, the deviations will also decrease, however at a much increased computational cost.

In total, the asymmetric model contains up to seven free parameters, namely R_p, h_cloud to describe the form of the planet, a_S, b, φ₀ to describe the orbit of the planet, and u₁, u₂ to describe the stellar intensity distribution.

4. Cloud diagnostics

4.1. Bisector span

The bisector span is the first potential diagnostic of an asymmetric transit that could indicate the presence of clouds. The bisector B as a function of transit depth d is a geometrical line. It is defined as the center of a horizontal line that intersects the transit light curve at a given d (see Fig. 8 for an illustration).

For a symmetric transit of a uniform planet, the bisector is a straight, vertical line, centered at zero. However, for asymmetric planets, two phenomena are observed. Firstly, the bisector is non-zero, shifted by a few minutes compared to the cloud-free case. Secondly, it is no longer a straight line, rather, with transit depth, a triangular-like shape can be seen. This shape is caused by the varying slope of the light curve as cloudy and cloud-free hemispheres, respectively, begin/end their passage in front of the star.

Figure 9 illustrates the bisector span for different values of cloud altitudes (zero cloud altitude corresponding to a symmetric transit), as calculated by the asymmetric model (Sect. 3). The adopted planetary and stellar parameters were those for Kepler-7b. We note the noise-like features at high transit depth, where the transit light curve is very flat (d ≳ 0.008 for the blue line). These features are due to the finite numerical resolution in the model.

The bisector is strongly affected by the observational noise at high transit depths, where the light curve is flat. Hence small variations due to noise will drastically change the bisector span. Figure 10 (upper panel) shows the same bisector spans, but with a synthetic noise of 34 ppm added to the light curve. This is the expected noise for stars of magnitude 8 with PlaTO 2.0 (Rauer et al. 2014). The bisector variations with transit depth are still clearly visible.

When adding a more realistic noise with an amplitude of 500 ppm, comparable to Kepler magnitude 11–12 stars, the shape of the bisector is hard to analyze. It may be possible to exclude high cloud altitudes of the order of 20 000 km, but note, however, that such an altitude would correspond to about 13 scale heights in the case of Kepler-7b.

Two main characteristics are proposed to analyze the bisector, as a function of cloud altitude. These are, firstly, the minimum bisector at half ingress (when the entire leading hemisphere already covers the star and the trailing hemisphere begins covering the star), B_min. Since an offset in Φ₀ of a symmetric transit also leads to a non-zero B_min, the second criterion evaluates the shape of the bisector. It uses the difference ΔB between the bisector at the beginning of the ingress and at half ingress. Figure 11 shows these two diagnostics as a function of cloud altitude for a Kepler-7 setup. It is clearly seen that there is an almost linear relationship between the diagnostics and cloud altitude. The slope of the relationship depends on stellar characteristics (limb darkening) and the planetary orbit (especially a_S).

To develop a detection criterion, the rms of the bisector σ_B is calculated as the standard deviation over ingress and egress. This excludes the full-transit phase where the planet covers the star completely and the bisector is most sensitive to observational noise. Then, the median B_med over the same range is calculated to establish a non-zero bisector: ${\mathit{\sigma }}_{\mathrm{0}}\mathrm{=}\frac{{\mathit{B}}_{\mathrm{med}}}{{\mathit{\sigma }}_{\mathrm{B}}}\mathrm{·}$(12)Once a non-zero bisector has been established, the shape is analyzed with the following criterion on a smoothed bisector: ${\mathit{\sigma }}_{\mathrm{shape}}\mathrm{=}\frac{\mathrm{\Delta }\mathit{B}}{{\mathit{\sigma }}_{\mathrm{B}}}\mathrm{·}$(13)For the 34 ppm cases in Fig. 10, we thus find significant detections (σ₀ between 2 and 6, and subsequently, σ_shape between 3 and 5). However, in the 500 ppm cases, all but the 20 000 km bisector span are compatible with zero. Even the most favorable case only shows a σ₀ of the order of 1.5, and when smoothing the bisector, one finds a σ_shape around 4. This would imply a somewhat significant detection, but as mentioned above, corresponding to much more than 10 scale heights.

4.2. Residuals

The second cloud diagnostics are the residuals of the data with respect to a symmetric transit model, such as the model of Mandel & Agol (2002). An asymmetric transit will imprint a clear structure on the residuals during ingress and egress when compared to a symmetric model. Figure 12 shows the residuals between symmetric and asymmetric models, again, as above, for Kepler-7 system parameters. We compare asymmetric models with R_P = 18 R_⊕ and varying cloud altitudes to symmetric models with a planetary radius R_eq such that the asymmetric planet and the symmetric planet cover the same surface area: ${\mathit{R}}_{\mathrm{eq}}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}{\mathit{R}}_{\mathrm{P}}^{\mathrm{2}}\mathrm{+}\mathrm{(}{\mathit{R}}_{\mathrm{P}}\mathrm{+}{\mathit{h}}_{\mathrm{cloud}}{\mathrm{)}}^{\mathrm{2}}\mathrm{)}\mathit{.}$(14)It is clearly seen that in the short ingress and egress phases of the light curve, the residuals show distinct structures. As above for the bisector span, these residuals are somewhat sensitive to the observational noise. As illustrated in Fig. 13, for a synthetic Gaussian noise of 34 ppm, the residuals still clearly show a structure. Upon increasing the noise, however, even at the adopted cloud altitude of 10 000 km, the residuals do not show a sign of clouds anymore.

4.3. Degeneracies

Equation (14) shows the principal correlation between the parameters of the asymmetric model: in the absence of limb darkening, asymmetric and symmetric models will be identical during transit, if Eq. (14) is satisfied. If limb darkening is present, small differences between various parameter combinations allows for resolving these degeneracies. However, observational noise will tend to counteract this effect. Therefore, in the presence of significant noise, strong correlation between parameters is expected.

Combinations of planetary radius and cloud altitude that satisfy Eq. (14), however, will show slight timing variations as to the start and end of the transit ingress and egress phases. These different starting/ending times of the planetary transit can be compensated for by choice of Φ₀, thus leading to the second possible correlation between model parameters if the light curve is significantly affected by noise.

In Fig. 14, several residuals are shown with respect to a symmetric model. Note in particular the smaller vertical scale compared to Figs. 12 and 13.

Adding synthetic white noise to the model light curves used in Fig. 14, we were able to determine the maximum noise level, up to which a clear distinction between models, hence decorrelating the parameters, would still be possible. In general, a noise level of lower than about 50–100 ppm allows to clearly retrieve the original model.

5. Inverse model

We use the Bayesian formalism to calculate posterior probability values p(V_P | D) for the parameter vector V_P in the model, given a set D of observations. $\mathit{p}\mathrm{\left(}{\mathit{V}}_{\mathrm{P}}\mathrm{\right|}\mathit{D}\mathrm{)}\mathrm{\propto }\mathit{p}\mathrm{(}\mathit{D}\mathrm{\left|}{\mathit{V}}_{\mathrm{P}}\mathrm{\right)}\mathrm{·}\mathit{p}\mathrm{\left(}{\mathit{V}}_{\mathrm{P}}\mathrm{\right)}\mathit{.}$(15)The likelihood p(D | V_P) is calculated assuming independent measurements and identically-distributed Gaussian errors for the individual data points. The priors p(V_P) are taken to be uninformative over the entire parameter range allowed.

5.1. MCMC algorithm

To sample the full parameter space, we adopt a Markov Chain Monte Carlo (MCMC) approach. We use the emcee python package developed by Foreman-Mackey et al. (2013), which implements an algorithm described in Goodman & Weare (2010). emcee uses multiple chains (in this work, typically around 500) to sample the parameter space. In each step, for each chain, the algorithm proposes a new positions based on the position of the entire ensemble of chains. Compared to more traditional MCMC approaches such as Metropolis-Hastings, emcee converges much quicker and is less dependent on initial conditions. Furthermore, emcee does not need a tuning of the parameter of the proposal function.

5.2. Fitting procedure

We proceeded in two different steps for the MCMC simulations.

In a first step, a low numerical resolution in the forward model was used that is however sufficient to constrain stellar and orbital characteristics (i.e., a_S, b, u₁, u₂) to a high degree of confidence. To ensure good convergence and avoid any contamination by initial conditions, the chains were run for 500 steps (>5 auto-correlation lengths for each parameter). The first 1–2 auto-correlation lengths were considered as burn-in and discarded for the calculation of parameter uncertainties. Convergence was checked by inspecting visually the evolution of the mean of the entire ensemble. Initial positions were obtained with a random sample within the assumed prior to allow the sampler to start by exploring the entire parameter space. Uncertainty ranges are calculated by marginalizing over the posterior distribution, thinned by the auto-correlation length of the particular parameter in question. We then determine 68% and 95% credibility regions as the [0.16, 0.84] and [0.03, 0.98] median-centered percentiles, respectively, of the cumulative probability distributions (CDF). If the parameter distribution were to be Gaussian, these credibility regions would correspond to the 1 and 2σ uncertainties, respectively.

Then, in a second step, the 95% credibility regions obtained in the first step were used to initialize a high-resolution forward model. The resolution is increased such that the stellar disk is resolved better than a planetary scale height, and thus, the model becomes sensitive to the effect of clouds on the transit shape. However, this increases the computational cost by about a factor of 100.

5.3. Scenarios

For each planet, we considered three scenarios. The first scenario (five parameters) is a symmetric scenario (no clouds) with Φ₀ fixed at 0. The second scenario relaxes this constraint and fits additionally for Φ₀, but retains the assumption of no clouds, hence symmetric transits. The third scenario is fully asymmetric, fitting for all seven parameters in the model.

6. Results

6.1. Bisector spans

In Figs. 15–17, the bisector spans of the three target planets are shown. As expected from the data quality in Figs. 4–6, the bisector spans are much better constrained for HAT-P-7b and HD 209458b.

We define the scale height time, T_H, as the characteristic time scale on which cloud phenomena will impact the transit light curve. Essentially, it describes the time needed for the planet to advance by one scale height while transiting: ${\mathit{T}}_{\mathrm{H}}\mathrm{=}\frac{\mathit{H}}{{\mathit{v}}_{\mathrm{orb}}}\mathit{,}$(16)with H the scale height and v_orb the orbital velocity. For Kepler-7b, we find T_H = 10.7 s, whereas for the relatively high-gravity planets HAT-P-7b and HD 209458b, T_H = 2.7 s and T_H = 4.5 s, respectively. In Figs. 15–17, we show the scale height times T_H for the different planets (green lines: ±T_H, blue lines ±5T_H) to illustrate the potential constraints that could be inferred from the bisector spans.

It is clear that neither Kepler-7b nor HAT-P-7b show a definitive sign of asymmetry in the bisector span. Following the detection criteria developed in Eqs. (12) and (13), we find σ₀ of far less than unity for both Kepler-7b and HAT-P-7b. This is because for Kepler-7b the photometric noise is too high and for HAT-P-7b, T_H simply is too small. It should be noted, though, that for a Kepler-7b-like planet orbiting around a star similar to HAT-P-7 or HD 209458, the bisector span could at least have been used to constrain clouds to within a few scale heights, since Kepler-7b has a very favorable scale height.

In the case of HD 209458b, the bisector span is very well constrained to within a few scale heights. However, due to the low time sampling of the transit light curve of HD 209458b (only four transits, covered by a few hundred points, Brown et al. 2001), the bisector span of HD 209458b is not conclusive to infer reliably the presence of any asymmetry.

6.2. Residuals

In Appendix A, the parameter space projections of the symmetric MCMC simulations are shown (Figs. A.1, A.2, A.4 and A.5). All stellar, orbital and planetary parameters are well constrained for both planets. In the case of Kepler-7b, the symmetric models show that the observations are consistent with Φ₀ = 0, whereas for HAT-P-7b, a statistically significant deviation of the order of 2 s is found by the MCMC model. This is slightly larger than the 1σ uncertainty associated with T₀ in Table 1. The best-fitting models have a ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}$ value of 1 and 0.99 for Kepler-7b and HAT-P-7b, respectively. Hence the quality of the fit is satisfactory.

In Figs. 18 and 19, the residuals of the symmetric fits to the data are shown. As before with the bisector spans, the residuals do not show any hint of asymmetry or clouds comparable to the ones discussed above in Sect. 4.

Again, when comparing the results for both planets, it is apparent that a Kepler-7b-like planet around a relatively bright star such as HAT-P-7 would be much more favorable for such an analysis and could be amenable to cloud altitude constraints with this technique.

6.3. Asymmetric transit models

The parameter space projections of the asymmetric MCMC simulations are shown in Figs. A.3 and A.6). It is evident that in the case of Kepler-7b, no clear conclusions can be drawn from the analysis, since planetary radius R_p, cloud altitude h_cloud and phase offset φ₀ are strongly correlated parameters, as expected from the high noise level (see discussion above). Hence, for Kepler-7b, the analysis is compatible with h_cloud = 0.

Note, however, that in the case of HAT-P-7b, where the standard rms of the final light curve is 150 ppm instead of the roughly 600 ppm for Kepler-7b, the degeneracy is somewhat weaker, and some constraints appear to be feasible. From the marginalized posterior distributions, it seems that a slight asymmetry is preferred, with h_cloud< 0 at 68% confidence (although h_cloud = 0 is within the 95% credibility region). This would suggest that the trailing limb appears slightly larger (of the order of a couple of scale heights). Such a result is actually consistent with the analysis of von Paris et al. (2016) who found that the asymmetric phase curve of HAT-P-7b is more likely due to an asymmetry in scattering properties than due to a thermal offset, contrary to previous studies (e.g., Esteves et al. 2015). However, given that the constraints are very weak, such a conclusion appears to be only tentative.

We note in addition that the asymmetric model contains more parameters than the symmetric models. Therefore, we apply the Bayesian Information Criterion (BIC) to choose between models. The HAT-P-7b asymmetric model is preferred by a factor of about 400 over the purely symmetric model. However, the χ² of both the asymmetric model and the free-Φ₀ model are virtually identical, so this model in turn is preferred by about a factor of 400 over the asymmetric model (ΔBIC of about 12, respectively). Table 2 summarizes the fit results for the two planets and shows the planetary parameters (orbital and stellar parameters are virtually identical for all fit scenarios).

7. Discussion

7.1. Suitable targets

Although it affects the entire transit light curve via limb darkening, most of the effect of the phenomenon that is looked for occurs on extremely short timescales, of the order of a few to a few tens of seconds. Therefore, the light curve must be obtained at a very high cadence. This can be done with short exposure times or with a long time series covering dozens, or even hundreds, of orbits. Furthermore, bright targets are needed to obtain the required photometric precision of the order of a few tens of ppm.

In addition to the requirements for the star, the planets must also be suitable for this technique to be effective. Orbital speeds for close-in planets are of the order of 100–200 km s^-1, hence the scale height of the considered planet should be of the order of 1000–2000 km. Low-gravity planets in short orbits have very large scale heights ($\mathit{H}\mathrm{~}\frac{{\mathit{T}}_{\mathrm{eq}}}{\mathit{g}}$, g gravity, T_eq equilibrium temperature). For scale heights of the order of 10³–10⁴ km, the signature of clouds during primary transit will be visible for bright targets.

Figure 20 shows a compilation of transiting planets discovered so far. The host star V magnitude is used as a color coding, and we also show transit depth (x-axis) and surface gravity (y-axis). As is clear, suitable targets (in terms of the diagram, blue dots in the lower quarter of the diagram) are missing. For instance, as already noted above, HAT-P-7 is a reasonably bright target, but the surface gravity is much too large. Inversely, Kepler-7b has a very suitable surface gravity (only slightly larger than Mars), but orbits a too faint host star for this technique to work.

Many transiting planets have been discovered to date with very low surface gravities, down to ~1 m s^-2 (e.g., Hartman et al. 2011; Lissauer et al. 2011; Bonfils et al. 2012; Jontof-Hutter et al. 2014). However, as for Kepler-7b, their host stars are too faint for cloud detection with primary transits.

Upcoming space missions such as ChEOPS (Broeg et al. 2013), TESS (Ricker et al. 2014) and PlaTO 2.0 (Rauer et al. 2014) are expected to deliver many more transiting planets and candidates as well as high-precision photometry of radial-velocity planets, but, in contrast to CoRoT and Kepler, around much brighter host stars (down to magnitude 8). Therefore, the detection of cloud signatures via primary transit asymmetries should be feasible for at least a few candidates from these surveys.

7.2. Forward scattering clouds

Optically thick clouds are responsible for the reflected light contribution to the observed phase curves and are the aim of this work. However, optically thin clouds can also play a large role in shaping the transmission across the planetary limbs. Depending on cloud particle characteristics, a strong forward scattering peak can appear (e.g., Wakeford & Sing 2015), which then scatters stellar light into the line-of-sight. Therefore, the presence of high-altitude, optically thin clouds can lead to a decreased contrast between cloud-free and cloud-covered limbs. Detailed modeling of such situations is however beyond the scope of this paper.

7.3. Other sources of asymmetric transits

Uniquely attributing asymmetric transit light curves to inhomogeneous cloud cover is not possible. Over the last decade, several other astrophysical sources of asymmetric transit light curves have been presented.

For instance, 3D atmospheric modeling of hot Jupiters generally predicts strong eastwards circulation (e.g., Showman & Guillot 2002), and thus, a hotter trailing limb than the morning limb (e.g., Agúndez et al. 2012, 2014; Parmentier et al. 2013). It follows that the atmosphere extends further, hence the transit depth is larger. Agúndez et al. (2014) note that this thermal effect dominates over the effect of chemistry in transmission spectra. In the science case considered in this work, the hotter evening side somewhat diminishes the effect of the inhomogeneous cloud cover, since the “missing” planetary radius is compensated for by the extended atmosphere on the trailing limb. On the other hand, for otherwise homogenous planets (for instance, with global cloud cover), the purely thermal effect could produce slightly asymmetric transits.

In an order-of-magnitude estimate, when calculating the extent of the atmosphere from the homogenized high-pressure layer at p_H to the cloud layer at p_C, one finds that the evening atmosphere is extended (compared to the morning atmosphere) by an amount h_E, ${\mathit{h}}_{\mathrm{E}}\mathrm{=}\ln \left(\frac{{\mathit{p}}_{\mathrm{H}}}{{\mathit{p}}_{\mathrm{C}}}\right)\frac{\mathrm{\Delta }\mathit{T}}{{\mathit{T}}_{\mathrm{M}}}\mathrm{·}\mathit{H,}$(17)where H is the scale height, ΔT is the temperature difference between morning temperature T_M and evening temperature T_E. For many planets, 3D circulation models find relatively moderate values of about $\frac{\mathrm{\Delta }\mathit{T}}{{\mathit{T}}_{\mathrm{M}}}\mathrm{\approx }\mathrm{10}$–20%. Using p_H = 1 bar and p_C = 50 mbar as representative values, one then finds h_E ≈ 0.3–0.6 H. This is not negligible, however not expected to be the dominant photometric effect.

Other possible physical causes of asymmetric transits are, e.g., the presence of exomoons around the planet in question (e.g., Szabó et al. 2006; Kipping 2009; Heller 2014) or gas and dust tails of strongly irradiated or disrupting planets (e.g., Brogi et al. 2012; Rappaport et al. 2014; Sanchis-Ojeda et al. 2015).

A possible way to discriminate between such scenarios would be, for instance, that in the case of tailed exoplanets, the transit depth would vary quite remarkably. Previous studies of HAT-P-7b have revealed a slight variation of the transit depth with the observation quarter (Van Eylen et al. 2013), which however is not attributed to physical reasons, but rather to instrumental effects. Detailed model comparisons between competing scenarios are however not the subject of this work, but rather for future studies.

8. Conclusions

We have presented a new technique for identifying cloud signatures in the transit light curves of exoplanets. Inhomogeneous cloud cover produces asymmetries in the transit light curve. Such asymmetries reveal themselves via the bisector span and characteristic residuals that differ from the standard, symmetric models.

We have applied this reasoning to search for hints of clouds in the transits of Kepler-7b, HAT-P-7b, and HD 209458b. No clear signs of clouds have been found, although for HAT-P-7b, results seem to indicate a small but not statistically very significant asymmetry consistent with previous phase curve analyses by von Paris et al. (2016). There are several reasons for this. In the case of Kepler-7b, the photometric noise is too high. HAT-P-7b, despite its better photometric quality, has a scale height which is too small to allow constraints on the presence of clouds or the cloud altitude. In the case of HD 209458b, the time sampling of the used transit light curve is too small, which again prevents a clear detection of clouds.

However, we have shown that with future photometric surveys and suitable targets, asymmetric transits will be detectable. This will help to better constrain cloud properties on exoplanets.

---

Acknowledgments

This study has received financial support from the French State in the frame of the “Investments for the future” Programme IdEx Bordeaux, reference ANR-10-IDEX-03-02. Computer time for this study was provided in part by the computing facilities MCIA (Mésocentre de Calcul Intensif Aquitain) of the Université de Bordeaux and of the Université de Pau et des Pays de l’Adour. P.G. acknowledges support from the ERC Starting Grant (3DICE, grant agreement 336474). P.G.’s postdoctoral position is funded by INSU/CNRS. This paper includes data collected by the Kepler mission. Funding for the Kepler mission is provided by the NASA Science Mission directorate. Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NNX09AF08G and by other grants and contracts. This research has made use of the Exoplanet Orbit Database and the Exoplanet Data Explorer at exoplanets.org.

References

Appendix A: MCMC results

All Tables

All Figures