# Inferring asymmetric limb cloudiness on exoplanets from transit light curves

P. von Paris, P. Gratier, J. Leconte, Franck Selsis, Pascal Bordé

► To cite this version:

P. von Paris, P. Gratier, J. Leconte, Franck Selsis, Pascal Bordé. Inferring asymmetric limb cloudiness on exoplanets from transit light curves. Astronomy & Astrophysics - A&A, 2016, 589, pp.id.A52. 10.1051/0004-6361/201527894. [hal-01275233](https://hal.science/hal-01275233v1)

HAL Id: hal-01275233

https://hal.science/hal-01275233v1

Submitted on 11 Nov 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Inferring asymmetric limb cloudiness on exoplanets from transit light curves

P. von Paris<fnref n="1" />,<fnref n="2" />, P. Gratier<fnref n="1" />,<fnref n="2" />, P. Bordé<fnref n="1" />,<fnref n="2" />, J. Leconte<fnref n="1" />,<fnref n="2" />, and F. Selsis<fnref n="1" />,<fnref n="2" />

Received 4 December 2015 / Accepted 12 February 2016

## ABSTRACT

*Context.* Clouds have been shown to be present in many exoplanetary atmospheres. Cloud formation modeling predicts considerable inhomogeneities of cloud cover, consistent with optical phase curve observations. However, optical phase curves cannot resolve some existing degeneracies between cloud location and cloud optical properties.

*Aims.* We present a conceptually simple technique for detecting inhomogeneous cloud cover on exoplanets. Such an inhomogeneous cloud cover produces an asymmetric primary transit of the planet in front of the host star. Asymmetric transits produce characteristic residuals that are different from standard symmetric models. Furthermore, bisector spans can be used to determine asymmetries in the transit light curve.

*Methods.* We apply a model of asymmetric transits to the light curves of HAT-P-7b, Kepler-7b, and HD 209458b and search for possible cloud signatures. The nearly uninterrupted *Kepler* photometry is particularly well suited for this method since it allows for a very high time resolution.

*Results.* We do not find any statistically sound cloud signature in the data of the considered planets. For HAT-P-7b, a tentative detection of an asymmetric cloud cover is found, consistent with analysis of the optical phase curve. Based on Bayesian probability arguments, a symmetric model with an offset in the transit ephemeris is still the most viable model. This work demonstrates that for suitable targets, namely low-gravity planets around bright stars, the method can be used to constrain cloud cover characteristics and is thus a helpful additional tool for the study of exoplanetary atmospheres.

**Key words.** techniques: photometric – planets and satellites: atmospheres – planets and satellites: individual: Kepler-7b – planets and satellites: individual: HAT-P-7b – planets and satellites: individual: HD 209458b

## 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 black-body 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

<footnote n="1">¹ Université Bordeaux, LAB, UMR 5804, 33270 Floirac, France
e-mail: vonparisp@gmail.com</footnote>

<footnote n="2">² CNRS, LAB, UMR 5804, 33270 Floirac, France</footnote>
“morning” side across the substellar meridian to the hotter “evening” side<fnref n="1" />.

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<fnref n="2" />. 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

<figure><img src="image_2.png" /><figcaption>Fig. 1. SC-SAP for Kepler-7 in quarter Q3.</figcaption></figure>

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<fnref n="3" />. 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,

$$M_{2.4}(j) = \text{median}(R(i), |t(i) - t(j)| \le 1.2 \text{ h}), \quad (1)$$

where $t$ denotes the time stamps of the data point.

From there, we calculate the median absolute deviation (MAD) as follows:

$$\text{MAD} = \text{median}(|R - M_{2.4}|). \quad (2)$$

A point $j$ is rejected if

$$|R(j) - M_{2.4}(j)| \ge 4 \times \text{MAD}. \quad (3)$$

Figure 2 illustrates this procedure, for the 9th transit of Kepler-7b in Q3 (green lines).

2. Normalization: After outlier removal, we calculate a running median $M_{48}$ for each point $j$, using a roughly 2-day window (red line in Fig. 2):

$$M_{48}(j) = \text{median}(R(i), |t(i) - t(j)| \le 24 \text{ h}). \quad (4)$$

The normalized light curve $N$ is obtained by dividing the raw light curve by the median-filtered light curve:

$$N = \frac{R}{M_{48}}. \quad (5)$$

<footnote n="1">¹ Morning and evening refer here to an assumed global eastward atmospheric circulation.</footnote>

<footnote n="2">² While submitting this manuscript for publication, we became aware of another paper that also advocates this idea (Line & Parmentier 2016). However, they focus on the transmission spectra and do not concentrate further on photometric data.</footnote>

<footnote n="3">³ http://archive.stsci.edu/kepler/data_search/search.php</footnote>
<figure><img src="image_3.png" /><figcaption>**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).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Normalized SC-SAP for Kepler-7 in Q3.</figcaption></figure>

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.

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,

$$ \Phi(j) = \frac{t(j) - T_0}{P} - \lfloor \frac{t(j) - T_0}{P} \rfloor, \qquad (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)<fnref n="4" /> and summarized in Table 1.

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Phase-folded SC-SAP for Kepler-7 in Q3-Q8.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Phase-folded SC-SAP for HAT-P-7b in Q0-Q17.</figcaption></figure>

<table><caption>**Table 1.** DR24 values for *T*₀ and *P*, used in Eq. (6).</caption><thead><tr><td>Planet</td><td><i>T</i><sub>0</sub>-2 454 833 [days]</td><td><i>P</i> [days]</td></tr></thead><tbody><tr><td>Kepler-7b</td><td>134.2768383</td><td>4.885488953</td></tr><tr><td>HAT-P-7b</td><td>121.3585723</td><td>2.204735365</td></tr></tbody></table>

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*<sub>trail</sub> and *R*<sub>lead</sub> for the trailing and leading

<footnote n="4">⁴ http://exoplanetarchive.ipac.caltech.edu</footnote>
<figure><img src="image_7.png" /><figcaption>Fig. 6. Transit light curve of HD 209458b. Data from Brown et al. (2001).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Illustration of an asymmetric transit: the leading, cloudy hemisphere has a larger radius than the trailing, cloud-free hemisphere.</figcaption></figure>

hemisphere, respectively. These two radii are parameterized as follows:

$$R_{\text{trail}} = R_p \qquad (7)$$

$$R_{\text{lead}} = R_p + h_{\text{cloud}}, \qquad (8)$$

where $R_p$ is the planetary radius and $h_{\text{cloud}}$ is the cloud altitude. We note that $h_{\text{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_* \times 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_* = 10000$ 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_1$ and $u_2$ as the linear and the quadratic limb-darkening coefficients for the intensity distribution $I$ across the stellar disk:

$$I(\mu) = 1 - u_1 \cdot (1 - \mu) - u_2 \cdot (1 - \mu)^2, \qquad (9)$$

with $\mu = \cos(\phi)$ and $\phi$ 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:

$$a_S = \frac{a}{R_*}, \qquad (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:

$$b = \cos(i) \cdot a_S. \qquad (11)$$

In addition, we introduce a seventh parameter which is the phase offset $\phi_0$ 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 $\phi_0 = 0$, however due to, e.g., eccentric orbits, $\phi_0$ 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_{\text{cloud}}$ to describe the form of the planet, $a_S$, $b$, $\phi_0$ to describe the orbit of the planet, and $u_1, u_2$ 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 \gtrsim 0.008$ for the blue line). These features are due to the finite numerical resolution in the model.
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Definition of the bisector *B* as a function of depth *d*. Symmetric transit in black, asymmetric transit with “morning” clouds in red.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Bisector values for different adopted $h_{cloud}$. No observational noise present.</figcaption></figure>

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 bisec-
tor, as a function of cloud altitude. These are, firstly, the mini-
mum bisector at half ingress (when the entire leading hemisphere
already covers the star and the trailing hemisphere begins cover-
ing the star), $B_{min}$. Since an offset in $\Phi_0$ 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 $\Delta B$ between the bisec-
tor at the beginning of the ingress and at half ingress. Figure 11
shows these two diagnostics as a function of cloud altitude for a

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Same as Fig. 9, but with synthetic Gaussian noise on the original light curve. *Upper panel:* 34 ppm. *Lower panel:* 500 ppm. Note the differences in horizontal scale.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Bisector diagnostics as a function of cloud altitude for the same setup as in Fig. 9 (see text for discussion).</figcaption></figure>

Kepler-7 setup. It is clearly seen that there is an almost linear re-
lationship between the diagnostics and cloud altitude. The slope
<figure><img src="image_14.png" /><figcaption>Fig. 12. Residuals between symmetric and asymmetric models.</figcaption></figure>

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 $\sigma_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:

$$ \sigma_0 = \frac{B_{med}}{\sigma_B}. \qquad (12) $$

Once a non-zero bisector has been established, the shape is analyzed with the following criterion on a smoothed bisector:

$$ \sigma_{\text{shape}} = \frac{\Delta B}{\sigma_B}. \qquad (13) $$

For the 34 ppm cases in Fig. 10, we thus find significant detections ($\sigma_0$ between 2 and 6, and subsequently, $\sigma_{\text{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 $\sigma_0$ of the order of 1.5, and when smoothing the bisector, one finds a $\sigma_{\text{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_{\text{eq}}$ such that the asymmetric planet and the symmetric planet cover the same surface area:

$$ R_{\text{eq}}^2 = \frac{1}{2} (R_P^2 + (R_P + h_{\text{cloud}})^2). \qquad (14) $$

It is clearly seen that in the short ingress and egress phases of the light curve, the residuals show distinct structures. As above

<figure><img src="image_15.png" /><figcaption>Fig. 13. Residuals between symmetric and asymmetric models with a cloud altitude of 10 000 km: noise levels of 34 and 500 ppm, respectively.</figcaption></figure>

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 $\Phi_0$, 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.

$$ p(V_P|D) \propto p(D|V_P) \cdot p(V_P). \qquad (15) $$
<figure><img src="image_16.png" /><figcaption>Fig. 14. Residuals between various asymmetric models and a symmetric model ($\Phi_0 = 0, h_{cloud} = 0$).</figcaption></figure>

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_1$, $u_2$) 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\sigma$ uncertainties, respectively.

<figure><img src="image_17.png" /><figcaption>Fig. 15. Kepler-7b bisector (Q3-Q8). Red horizontal line indicates zero. Green and blue dashed lines represent $\pm 1$ and $\pm 5T_H$, respectively (see text).</figcaption></figure>

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 $\Phi_0$ fixed at 0. The second scenario relaxes this constraint and fits additionally for $\Phi_0$, 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:

$$T_H = \frac{H}{v_{orb}}, \qquad (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: $\pm T_H$, blue lines $\pm 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 $\sigma_0$ of far
<figure><img src="image_18.png" /><figcaption>**Fig. 16.** HAT-P-7b bisector (Q0–Q17). Red horizontal line indicates zero. Green and blue dashed lines represent ±1 and ±5*T<sub>H</sub>*, respectively. Note the smaller horizontal scale compared to Fig. 15.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Fig. 17.** HD 209458b bisector. Red horizontal line indicates zero. Green and blue dashed lines represent ±1 and ±5*T<sub>H</sub>*, respectively. Note the smaller horizontal scale compared to Fig. 16.</figcaption></figure>

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*<sub>H</sub> 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

<figure><img src="image_20.png" /><figcaption>**Fig. 18.** Best-fit residuals for Kepler-7b.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Fig. 19.** Best-fit residuals for HAT-P-7b. Note the difference in vertical scale compared to Fig. 18.</figcaption></figure>

for both planets. In the case of Kepler-7b, the symmetric models show that the observations are consistent with $\Phi_0 = 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\sigma$ uncertainty associated with $T_0$ in Table 1. The best-fitting models have a $\chi^2_{red}$ 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
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 homogeneous 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$,

$$h_E = \ln \left( \frac{p_H}{p_C} \right) \frac{\Delta T}{T_M} \cdot H, \quad (17)$$

where $H$ is the scale height, $\Delta 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{\Delta T}{T_M} \approx 10-20\%$. Using $p_H = 1$ bar and $p_C = 50$ mbar as representative values, one then finds $h_E \approx 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.

*Acknowledgements.* 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

Agúndez, M., Venot, O., Iro, N., et al. 2012, A&A, 548, A73

Agúndez, M., Parmentier, V., Venot, O., Hersant, F., & Selsis, F. 2014, A&A, 564, A73

Barclay, T., Huber, D., Rowe, J. F., et al. 2012, ApJ, 761, 53

Bean, J. L., Kempton, E., & Homeier, D. 2010, Nature, 468, 669

Bonfils, X., Gillon, M., Udry, S., et al. 2012, A&A, 546, A27

Borucki, W. J., Koch, D., Jenkins, J., et al. 2009, Science, 325, 709

Broeg, C., Fortier, A., Ehrenreich, D., et al. 2013, in Eur. Phys. J. Web Conf., 47, 3005

Brogi, M., Keller, C. U., de Juan Ovelar, M., et al. 2012, A&A, 545, L5

Brown, T. M., Charbonneau, D., Gilliland, R. L., Noyes, R. W., & Burrows, A. 2001, ApJ, 552, 699

Deming, D., Wilkins, A., McCullough, P., et al. 2013, ApJ, 774, 95
Demory, B.-O., Seager, S., Madhusudhan, N., et al. 2011, ApJ, 735, L12

Demory, B.-O., de Wit, J., Lewis, N., et al. 2013, ApJ, 776, L25

Esteves, L. J., De Mooij, E. J. W., & Jayawardhana, R. 2013, ApJ, 772, 51

Esteves, L. J., De Mooij, E. J. W., & Jayawardhana, R. 2015, ApJ, 804, 150

Evans, T. M., Pont, F., Sing, D. K., et al. 2013, ApJ, 772, L16

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Gibson, N. P., Aigrain, S., Barstow, J. K., et al. 2013, MNRAS, 436, 2974

Goodman, J., & Weare, J. 2010, Commun. Appl. Math. Comput. Sci., 5

Han, E., Wang, S. X., Wright, J. T., et al. 2014, PASP, 126, 827

Hansen, C. J., Schwartz, J. C., & Cowan, N. B. 2014, MNRAS, 444, 3632

Hartman, J. D., Bakos, G. Á., Kipping, D. M., et al. 2011, ApJ, 728, 138

Heller, R. 2014, ApJ, 787, 14

Jontof-Hutter, D., Lissauer, J. J., Rowe, J. F., & Fabrycky, D. C. 2014, ApJ, 785, 15

Kipping, D. M. 2009, MNRAS, 392, 181

Knutson, H. A., Benneke, B., Deming, D., & Homeier, D. 2014, Nature, 505, 66

Kreidberg, L., Bean, J. L., Désert, J.-M., et al. 2014, Nature, 505, 69

Latham, D. W., Borucki, W. J., Koch, D. G., et al. 2010, ApJ, 713, L140

Lee, G., Helling, C., Dobbs-Dixon, I., & Juncher, D. 2015, A&A, 580, A12

Line, M. R., & Parmentier, V. 2016, ApJ, 820, 78

Lissauer, J. J., Fabrycky, D. C., Ford, E. B., et al. 2011, Nature, 470, 53

Mandel, K., & Agol, E. 2002, ApJ, 580, L171

Mazeh, T., & Faigler, S. 2010, A&A, 521, L59

Pál, A., Bakos, G. Á., Torres, G., et al. 2008, ApJ, 680, 1450

Parmentier, V., Showman, A. P., & Lian, Y. 2013, A&A, 558, A91

Quintana, E. V., Rowe, J. F., Barclay, T., et al. 2013, ApJ, 767, 137

Rappaport, S., Barclay, T., DeVore, J., et al. 2014, ApJ, 784, 40

Rauer, H., Catala, C., Aerts, C., et al. 2014, Exp. Astron., 38, 249

Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2014, SPIE Conf. Ser., 9143, 20

Sanchis-Ojeda, R., Rappaport, S., Pallè, E., et al. 2015, ApJ, 812, 112

Schlawin, E., Zhao, M., Teske, J. K., & Herter, T. 2014, ApJ, 783, 5

Showman, A. P., & Guillot, T. 2002, A&A, 385, 166

Sing, D. K., Wakeford, H. R., Showman, A. P., et al. 2015, MNRAS, 446, 2428

Snellen, I. A. G., de Mooij, E. J. W., & Albrecht, S. 2009, Nature, 459, 543

Sudarsky, D., Burrows, A., & Pinto, P. 2000, ApJ, 538, 885

Szabó, G. M., Szatmáry, K., Divéki, Z., & Simon, A. 2006, A&A, 450, 395

Van Eylen, V., Kjeldsen, H., Christensen-Dalsgaard, J., & Aerts, C. 2012, Astron. Nachr., 333, 1088

Van Eylen, V., Lindholm Nielsen, M., Hinrup, B., Tingley, B., & Kjeldsen, H. 2013, ApJ, 774, L19

von Paris, P., Gratier, P., Bordé, P., & Selsis, F. 2016, A&A, 587, A149

Wakeford, H. R., & Sing, D. K. 2015, A&A, 573, A122

Webber, M. W., Lewis, N. K., Marley, M., et al. 2015, ApJ, 804, 94

Welsh, W. F., Orosz, J. A., Seager, S., et al. 2010, ApJ, 713, L145

Winn, J. N., Johnson, J. A., Albrecht, S., et al. 2009, ApJ, 703, L99

Wright, J. T., Fakhouri, O., Marcy, G. W., et al. 2011, PASP, 123, 412
Appendix A: MCMC results

<figure><img src="image_22.png" /><figcaption>Fig. A.1. Kepler-7b posterior projections for the symmetric model ($\Phi_0 = 0, h_{\text{cloud}} = 0$). Dashed vertical lines represent marginalized 95% credibility regions. Smoothed 68% and 95% credibility regions in dark and light grey shade, respectively.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**Fig. A.2.** Kepler-7b posterior projections for the symmetric model with free $\Phi_0$ and $h_{\text{cloud}} = 0$. Dashed vertical lines represent marginalized 95% credibility regions. Smoothed 68% and 95% credibility regions in dark and light grey shade, respectively.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>**Fig. A.3.** Kepler-7b posterior projections for the asymmetric model. Dashed vertical lines represent marginalized 95% credibility regions. Smoothed 68% and 95% credibility regions in dark and light grey shade, respectively.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>**Fig. A.4.** HAT-P-7b posterior projections for the symmetric model ($\Phi_0 = 0, h_{\text{cloud}} = 0$). Dashed vertical lines represent marginalized 95% credibility regions. Smoothed 68% and 95% credibility regions in dark and light grey shade, respectively.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>**Fig. A.5.** HAT-P-7b posterior projections for the symmetric model with free $\Phi_0$ and $h_{\text{cloud}} = 0$. Dashed vertical lines represent marginalized 95% credibility regions. Smoothed 68% and 95% credibility regions in dark and light grey shade, respectively.</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>**Fig. A.6.** HAT-P-7b posterior projections for the asymmetric model. Dashed vertical lines represent marginalized 95% credibility regions. Smoothed 68% and 95% credibility regions in dark and light grey shade, respectively.</figcaption></figure>