\section{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., \citealp{bean2010}, \citealp{kreidberg2014}), GJ436b (e.g., \citealp{knutson2014}), CoRoT-1b (e.g., \citealp{schlawin2014}), WASP-31b (e.g., \citealp{sing2015}) or HAT-P-32b (e.g., \citealp{gibson2013}), 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, which mask the expected molecular and atomic absorption.

During secondary eclipse, the apparent dayside brightness temperature can be measured at IR wavelengths. Again, as different wavelengths probe different atmospheric pressure due to wavelength-dependent opacity, 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., \citealp{hansen2014}). A single blackbody brightness temperature can be assigned to the photosphere, independent of wavelength. Again, 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., \citealp{evans2013}). 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., \citealp{snellen2009}, \citealp{mazeh2010}, \citealp{barclay2012}, \citealp{quintana2013}, \citealp{esteves2013,esteves2015}), 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., \citealp{demory2013inhomogen}, \citealp{esteves2015}, \citealp{webber2015}). These are conceptually interpreted to be 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\footnote{Morning and evening refer here to an assumed global eastward atmospheric circulation.}.

In recent years, detailed theoretical cloud modeling has been performed (e.g., \citealp{sudarsky2000}, \citealp{parmentier2013}, \citealp{wakeford2015}, \citealp{webber2015}, \citealp{lee2015}). 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\footnote{While submitting this manuscript for publication, we became aware of another paper that also advocates this idea \citep{line2015}. However, they focus on the transmission spectra and do not concentrate further on photometric data.}. 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., \citealp{latham2010}, \citealp{demory2011kepler7,demory2013inhomogen}, \citealp{esteves2015}) and HAT-P-7b (\citealp{pal2008}, \citealp{borucki2009}, \citealp{welsh2010}, \citealp{esteves2015}, \citealp{vparis2015phase}). Furthermore, we analyze transit light curves of HD209458b taken from \citet{brown2001} (their Figure 3) obtained from four transits observed with the Hubble Space Telescope. HD209458 is one of the brightest stars around which a transiting planet has been found, and the transmission spectrum of HD209458b shows some evidence of possible cloud and/or haze cover (e.g., \citealp{deming2013}). 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: Section \ref{datared} presents the data reduction process adapted for this work, and Sect. \ref{formodel} presents the numerical model used in the transit modeling. Section \ref{clouddiag} 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. \ref{invmodel}. Results are presented in Sect. \ref{results} and discussed in Sect. \ref{discuss}. We conclude with Sect. \ref{summary}.

\section{Data reduction}

\label{datared}

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\footnote{http://archive.stsci.edu/kepler/data\_search/search.php}. An example is shown in Fig. \ref{sap_q3}.

\begin{figure}[h]
\begin{center}
\includegraphics[width=250pt]{K7_Q3_sap_sc}\\
\caption{SC-SAP for Kepler-7 in quarter Q3.}
\label{sap_q3}
\end{center}
\end{figure}

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:

\begin{enumerate}

\item Outlier removal: We calculate a running median $M_{2.4}$ for each point $j$ of the light curve, using a 2.4-hour window.

\begin{equation}
\label{median1hour}
M_{2.4}(j)=\rm{median}\left(R(i), |t(i)-t(j)|\leqslant 1.2 \rm{hrs}\right),
\end{equation}

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

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

\begin{equation}
\label{mad_def}
\rm{MAD}=\rm{median}\left(|R-M_{2.4}|\right).
\end{equation}

A point $j$ is rejected if

\begin{equation}
\label{reject}
|R(j)-M_{2.4}(j)|\geqslant 4\times \rm{MAD}.
\end{equation}

Figure \ref{transit9} illustrates this procedure, for the 9th transit of Kepler-7b in Q3 (green lines).

\begin{figure}[h]
\begin{center}
\includegraphics[width=250pt]{K7_Q3_T9}\\
\caption{9th transit in Q3: Illustration of running 2-day median (red line) and 2.4-hour median (green plain line). Green dashed lines encompass the 4x$\rm{MAD}$ criterion to identify outliers (see text).}
\label{transit9}
\end{center}
\end{figure}

\item 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. \ref{transit9}).

\begin{equation}
\label{median1day}
M_{48}(j)=\rm{median}\left(R(i), |t(i)-t(j)|\leqslant 24 \rm{hrs}\right).
\end{equation}

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

\begin{equation}
\label{normalized}
N=\frac{R}{M_{48}}.
\end{equation}

In this way, slow variations (due to, e.g., planetary phase curve, stellar activity, etc.) are mostly removed from the final light curve. Figure \ref{normalq3} shows the normalized light curve for Kepler-7 in quarter 3.

\begin{figure}[h]
\begin{center}
\includegraphics[width=250pt]{K7_Q3_normalized}\\
\caption{Normalized SC-SAP for Kepler-7 in Q3.}
\label{normalq3}
\end{center}
\end{figure}

\item Phase folding: The normalized light curve $N$ is then phase-folded to obtain the final light curve $L$ as a function of orbital phase $\Phi$. For this, $\Phi$ is defined such that $\Phi$=0 at primary transit.

\begin{equation}
\label{phasedef}
\Phi(j)=\frac{t(j)-T_0}{P}-\lfloor \frac{t(j)-T_0}{P} \rfloor,
\end{equation}

where $T_0$ is the mid-transit time of the first transit, $P$ the orbital period and $\lfloor x \rfloor$ represents the floor function, i.e., the  greatest integer less than or equal to $x$. Figure \ref{phasefoldedq3} shows the final light curve of Kepler-7b, using all quarters with SC data.

\end{enumerate}

\begin{figure}[h]
\begin{center}
\includegraphics[width=250pt]{K7_Q3Q8_data}\\
\caption{Phase-folded SC-SAP for Kepler-7 in Q3-Q8.}
\label{phasefoldedq3}
\end{center}
\end{figure}

The values for $T_0$ and $P$ are taken from the latest data release (DR24)\footnote{http://exoplanetarchive.ipac.caltech.edu} and summarized in Table \ref{dr24}.

\begin{table}[h]
  \centering
  \caption{DR24 values for $T_0$ and $P$, used in eq. \ref{phasedef}.}\label{dr24}
  \begin{tabular}{l|ll}
\hline
\hline
   Planet& $T_0$-2454833 [days] & $P$ [days] \\
\hline
  Kepler-7b &134.2768383 &4.885488953\\
   HAT-P-7b & 121.3585723 & 2.204735365
\end{tabular}
\end{table}

Figure \ref{phasefoldedq3_hat} shows the Q0-Q17 light curve of HAT-P-7b. Since HAT-P-7 is about 2.5 magnitudes brighter than Kepler-7 in the Kepler bandpass, the resulting light curve is much cleaner.

\begin{figure}[h]
\begin{center}
\includegraphics[width=250pt]{H7_Q0Q17_data}\\
\caption{Phase-folded SC-SAP for HAT-P-7b in Q0-Q17.}
\label{phasefoldedq3_hat}
\end{center}
\end{figure}

Figure \ref{trans_hd} shows the transit light curve of HD209458b, used in this work. Data is taken from \citet{brown2001} (their Figure 3). Observations were made with the Hubble Space Telescope during four transits of HD209458b and then phase-folded to obtain a composite light curve. As is clearly seen, the photometric quality of the light curve is exceptional, since HD209458 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.

\begin{figure}[h]
\begin{center}
\includegraphics[width=250pt]{hd209485_data}\\
\caption{Transit light curve of HD209458b. Data from \citet{brown2001}.}
\label{trans_hd}
\end{center}
\end{figure}

\section{Forward model}

\label{formodel}

\label{toymodel}

In this model, the planet is assumed to be composed of two hemispheres with radius $R_{\rm{trail}}$ and $R_{\rm{lead}}$ for the trailing and leading hemisphere, respectively. These two radii are parameterized as follows:

\begin{eqnarray}\label{hemispheres}
 R_{\rm{trail}}& = & R_p \\
R_{\rm{lead}} & = & R_p+ h_{\rm{cloud}}
\end{eqnarray}

where $R_p$ is the planetary radius and $h_{\rm{cloud}}$ is the cloud altitude. Note that $h_{\rm{cloud}}$<0 is permitted, which then implies a larger trailing hemisphere. The concept of the model is illustrated in Fig. \ref{asymsketch}.

\begin{figure}[h]
\begin{center}
\includegraphics[width=250pt]{star}\\
\caption{Illustration of an asymmetric transit: The leading, cloudy hemisphere has a larger radius than the trailing, cloud-free hemisphere.}
\label{asymsketch}
\end{center}
\end{figure}

To calculate the stellar total luminosity and estimate the area obscured by the planet, we use a $N_{\ast}\times N_{\ast}$ grid to resolve the stellar disk, with its center corresponding to the stellar center (see Fig. \ref{asymsketch}). $N_{\ast}$ is chosen such that a stellar pixel corresponds to roughly half a scale height of the planet in question ($N_{\ast}$=4,000 for Kepler-7b, $N_{\ast}$=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_1$ and $u_2$ as the linear and the quadratic limb-darkening coefficients for the intensity distribution $I$ across the stellar disk:

\begin{equation}
\label{limbquad}
I(\mu)=1-u_1\cdot (1-\mu)- u_2\cdot (1-\mu)^2,
\end{equation}

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., \citealp{winn2009}, \citealp{demory2011kepler7,demory2013inhomogen}, \citealp{esteves2015}).

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:

\begin{equation}
\label{kdef}
a_S=\frac{a}{R_{\ast}},
\end{equation}

where $a$ is the semi-major axis and  $R_{\ast}$ is the stellar radius (fixed in this work due to strong asteroseismology constraints, e.g. \citealp{demory2011kepler7}, \citealp{eylen2012}).

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:

\begin{equation}
\label{impact_para}
b=\cos \left(i\right) \cdot a_S .
\end{equation}

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. \ref{asymsketch}) 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 \citet{mandel2002} 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. \ref{residu_illu}). 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_{\rm{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.

\section{Cloud diagnostics}

\label{clouddiag}

\subsection{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. \ref{bisect_def} for an illustration).

\begin{figure}[h]
\begin{center}
\includegraphics[width=250pt]{bisect_illu2}\\
\caption{Illustration of the definition of the bisector $B$ as a function of depth $d$. Symmetric transit in black, asymmetric transit with "morning" clouds in red.}
\label{bisect_def}
\end{center}
\end{figure}

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 \ref{bisect} 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. \ref{toymodel}). The adopted planetary and stellar parameters were those for Kepler-7b. 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.

\begin{figure}[h]
\begin{center}
\includegraphics[width=250pt]{bisector}\\
\caption{Bisector values for different adopted $h_{\rm{cloud}}$. No observational noise present.}
\label{bisect}
\end{center}
\end{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 \ref{bisect_noise} (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 \citep{rauer2014}. The bisector variations with transit depth are still clearly visible.

\begin{figure}[h]
\begin{center}
\includegraphics[width=250pt]{bisector_noise_34ppm}\\
\includegraphics[width=250pt]{bisector_noise_500ppm}\\
\caption{Same as Fig. \ref{bisect}, but with synthetic Gaussian noise on the original light curve. Upper panel: 34\,ppm. Lower panel: 500\,ppm. Note the differences in horizontal scale.}
\label{bisect_noise}
\end{center}
\end{figure}

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.

\begin{figure}[h]
\begin{center}
\includegraphics[width=250pt]{bisect_diag}\\
\caption{Bisector diagnostics as a function of cloud altitude for the same setup as in Fig. \ref{bisect} (see text for discussion).}
\label{bisect_diag}
\end{center}
\end{figure}

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_{\rm{min}}$. Since an offset in $\Phi_0$ of a symmetric transit also leads to a non-zero $B_{\rm{min}}$, the second criterion evaluates the shape of the bisector. It uses the difference $\Delta B$ between the bisector at the beginning of the ingress and at half ingress. Figure \ref{bisect_diag} 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 $\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_{\rm{med}}$ over the same range is calculated to establish a non-zero bisector:

\begin{equation}
\label{nonzerodetect}
\sigma_0=\frac{B_{\rm{med}}}{\sigma_B}.
\end{equation}

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

\begin{equation}
\label{shapedetect}
\sigma_{\rm{shape}}=\frac{\Delta B}{\sigma_B}.
\end{equation}

For the 34\,ppm cases in Fig. \ref{bisect_noise}, we thus find significant detections ($\sigma_0$ between 2 and 6, and subsequently, $\sigma_{\rm{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_{\rm{shape}}$ around 4. This would imply a somewhat significant detection, but as mentioned above, corresponding to much more than 10 scale heights.

\subsection{Residuals}

The second cloud diagnostics are the residuals of the data with respect to a symmetric transit model, such as the model of \citet{mandel2002}. An asymmetric transit will imprint a clear structure on the residuals during ingress and egress when compared to a symmetric model. Figure \ref{residu_illu} 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_{\oplus}$ and varying cloud altitudes to symmetric models with a planetary radius $R_{\rm{eq}}$ such that the asymmetric planet and the symmetric planet cover the same surface area:

\begin{equation}
\label{req_def}
R_{\rm{eq}}^2=\frac{1}{2}\left(R_P^2+(R_P+h_{\rm{cloud}})^2\right).
\end{equation}

\begin{figure}[h]
\begin{center}
\includegraphics[width=250pt]{illustrate_residual_eq}\\
\caption{Illustration of residuals between symmetric and asymmetric models.}
\label{residu_illu}
\end{center}
\end{figure}

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. \ref{residu_illu_noise}, 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.

\begin{figure}[h]
\begin{center}
\includegraphics[width=250pt]{illustrate_residual_eq_noise}\\
\caption{Illustration of residuals between symmetric and asymmetric models with a cloud altitude of 10,000\,km: noise levels of 34 and 500\,ppm, respectively.}
\label{residu_illu_noise}
\end{center}
\end{figure}

\subsection{Degeneracies}

Equation \ref{req_def} 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. \ref{req_def} 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. \ref{req_def}, 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. \ref{degen_illu}, several residuals are shown with respect to a symmetric model. Note in particular the smaller vertical scale compared to Figs. \ref{residu_illu} and \ref{residu_illu_noise}.

\begin{figure}[h]
\begin{center}
\includegraphics[width=250pt]{illustrate_degen}\\
\caption{Illustration of residuals between various asymmetric models and a symmetric model ($\Phi_0$=0, $h_{\rm{cloud}}$=0).}
\label{degen_illu}
\end{center}
\end{figure}

Adding synthetic white noise to the model light curves used in Fig. \ref{degen_illu}, 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.

\section{Inverse model}

\label{invmodel}

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.

\begin{equation}
\label{bayes}
p(V_P | D)\propto p(D | V_P) \cdot p(V_P).
\end{equation}

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.

\subsection{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 \citet{foreman2013}, which implements an algorithm described in \citet{goodman2010}. 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.

\subsection{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.

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.

\subsection{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.

\section{Results}

\label{results}

\subsection{Bisector spans}

In Figs. \ref{bisect_k7}-\ref{bisect_hd}, the bisector spans of the three target planets are shown. As expected from the data quality in Figs. \ref{phasefoldedq3}, \ref{phasefoldedq3_hat} and \ref{trans_hd}, the bisector spans are much better constrained for HAT-P-7b and HD209458b.

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:

\begin{equation}
\label{scaletime}
T_H=\frac{H}{v_{\rm{orb}}},
\end{equation}

with $H$ the scale height and $v_{\rm{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 HD209458b, $T_H$=2.7\,s and $T_H$=4.5\,s, respectively. In Figs. \ref{bisect_k7}-\ref{bisect_hd}, we show the scale height times $T_H$ for the different planets (green lines: $\pm$$T_H$, blue lines $\pm$5$T_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. \ref{nonzerodetect} and \ref{shapedetect}, we find $\sigma_{0}$ 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 HD209458, 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.

\begin{figure}[h]
\begin{center}
\includegraphics[width=250pt]{bisector_k7_q3q8}\\
\caption{Kepler-7b bisector (Q3-Q8). Red horizontal line indicates zero. Green and blue dashed lines represent $\pm$1 and $\pm$5$T_H$, respectively (see text).}
\label{bisect_k7}
\end{center}
\end{figure}

\begin{figure}[h]
\begin{center}
\includegraphics[width=250pt]{bisector_h7_q0q17}\\
\caption{HAT-P-7b bisector (Q0-Q17). Red horizontal line indicates zero. Green and blue dashed lines represent $\pm$1 and $\pm$5$T_H$, respectively. Note the smaller horizontal scale compared to Fig. \ref{bisect_k7}.}
\label{bisect_h7}
\end{center}
\end{figure}

In the case of HD209458b, 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 HD209458b (only four transits, covered by a few hundred points, \citealp{brown2001}), the bisector span of HD209458b is not conclusive to infer reliably the presence of any asymmetry.

\begin{figure}[h]
\begin{center}
\includegraphics[width=250pt]{bisector_hd209458}\\
\caption{HD209458b bisector. Red horizontal line indicates zero. Green and blue dashed lines represent $\pm$1 and $\pm$5$T_H$, respectively. Note the smaller horizontal scale compared to Fig. \ref{bisect_h7}.}
\label{bisect_hd}
\end{center}
\end{figure}

\subsection{Residuals}

In Appendix \ref{mcmcresults}, the parameter space projections of the symmetric MCMC simulations are shown (Figs. \ref{K7_triangle_no}, \ref{K7_triangle_t},  \ref{H7_triangle_no} and \ref{H7_triangle_t}). 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 $\Phi_0$=0, whereas for HAT-P-7b, a statistically significant deviation of the order of 2\,sec is found by the MCMC model. This is slightly larger than the 1\,$\sigma$ uncertainty associated with $T_0$ in Table \ref{dr24}. The best-fitting models have a $\chi^2_{\rm{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. \ref{showres} and \ref{showres2}, 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. \ref{clouddiag}.

\begin{figure}[h]
\begin{center}
\includegraphics[width=250pt]{K7_Q3Q8_transit_asym_dr24_high_res_no_h_no_t_residuals}\\
\caption{Best-fit residuals for Kepler-7b.}
\label{showres}
\end{center}
\end{figure}

\begin{figure}[h]
\begin{center}
  \includegraphics[width=250pt]{H7_Q0Q17_transit_asym_dr24_high_res_no_h_no_t_residuals}\\
\caption{Best-fit residuals for HAT-P-7b. Note the difference in vertical scale compared to Fig. \ref{showres}.}
\label{showres2}
\end{center}
\end{figure}

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.

\subsection{Asymmetric transit models}

The parameter space projections of the asymmetric MCMC simulations are shown in Figs. \ref{K7_triangle_asym} and \ref{H7_triangle_asym}).  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_{\rm{cloud}}$ and phase offset $\phi_0$ are strongly correlated parameters, as expected from the high noise level (see discussion above). Hence, for Kepler-7b, the analysis is compatible with $h_{\rm{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_{\rm{cloud}}$<0 at 68\,\% confidence (although $h_{\rm{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 \citet{vparis2015phase} 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., \citealp{esteves2015}). However, given that the constraints are very weak, such a conclusion appears to be only tentative.

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 $\chi^2$ of both the asymmetric model and the free-$\Phi_0$ model are virtually identical, so this model in turn is preferred by about a factor of 400 over the asymmetric model ($\Delta$BIC of about 12, respectively). Table \ref{mcmcsum} summarizes the fit results for the two planets and shows the planetary parameters (orbital and stellar parameters are virtually identical for all fit scenarios).

\begin{table*}
  \centering
  \caption{95\,\% credibility regions of the planetary parameters ($R_P$, $h_{\rm{cloud}}$, $\Phi_0$) of our model fits. $\Delta$BIC and model probability ratios $p_M$ with respect to best model are stated (best model in bold).}\label{mcmcsum}
  \begin{tabular}{l|llcc|lll}
\hline
\hline
   planet& scenario & $R_P$ [$R_J$]& $\Phi_0$ [min]& $h_{\rm{cloud}}$ [km] &  $\Delta$ BIC & $p_M$\\
\hline
\hline
Kepler-7b  &   &&&\\
  & \textbf{symmetric}  & 17.92<$R_P$<18.04&-&-&0&1\\
  & symmetric + $\Phi_0$  &17.92<$R_P$<18.04&-0.12<$\Phi_0$<0.05&-&-10.2&$6.0 \cdot 10^{-3}$\\
  & asymmetric  &16.74<$R_P$<18.91&-1.17<$\Phi_0$<0.90& -12,450<$h_{\rm{cloud}}$<15,370&-20.9&$2.8 \cdot 10^{-5}$\\
\hline
\hline
HAT-P-7b  &  && &&&\\
  & symmetric  &16.09<$R_P$<16.10&-&-&-23.9&$6.3 \cdot 10^{-6}$\\
 & \textbf{symmetric + $\Phi_0$}  &16.09<$R_P$<16.10&-0.03<$\Phi_0$<-0.01&- & 0&1\\
  & asymmetric  &15.98<$R_P$<16.39&-0.10<$\Phi_0$<0.19&-3,860<$h_{\rm{cloud}}$<1,370&-12.0& $2.4 \cdot 10^{-3}$\\
\end{tabular}
\end{table*}

\section{Discussion}

\label{discuss}

\subsection{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\,kms$^{-1}$, hence the scale height of the considered planet should be of the order of 1,000-2,000\,km. Low-gravity planets in short orbits have very large scale heights ($H\sim \frac{T_{\rm{eq}}}{g}$, $g$ gravity, $T_{\rm{eq}}$ equilibrium temperature). For scale heights of the order of 10$^3$-10$^4$\,km, the signature of clouds during primary transit will be visible for bright targets.

Figure \ref{targets_choice} 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 $\sim$1\,ms$^{-2}$ (e.g., \citealp{hartman2011}, \citealp{lissauer2011}, \citealp{bonfils2012}, \citealp{jontof2014}). However, as for Kepler-7b, their host stars are too faint for cloud detection with primary transits.

\begin{figure}[h]
\begin{center}
\includegraphics[width=250pt]{plot}\\
\caption{Exoplanet statistics: Transit depth (x-axis) and surface gravity ($\log g$, in cm s$^{-2}$, y-axis). Color coding for the host star V magnitude. Planets used in this work are indicated. Data and plot taken from exoplanets.org (see \citealp{wright2011} and \citealp{han2014} for detailed description of the database).}
\label{targets_choice}
\end{center}
\end{figure}

Upcoming space missions such as ChEOPS \citep{broeg2013}, TESS \citep{ricker2014} and PlaTO 2.0 \citep{rauer2014} 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.

\subsection{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., \citealp{wakeford2015}), 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.

\subsection{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., \citealp{showman2002}), and thus, a hotter trailing limb than the morning limb (e.g., \citealp{agundez2012,agundez2014}, \citealp{parmentier2013}). It follows that the atmosphere extends further, hence the transit depth is larger.  \citet{agundez2014} 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$:

\begin{equation}
\label{evening}
h_E=\ln \left(\frac{p_H}{p_C}\right) \frac{\Delta T}{T_M}\cdot H,
\end{equation}

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$=\,50mbar 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., \citealp{szabo2006}, \citealp{kipping2009ttv}, \citealp{heller2014}) or gas and dust tails of strongly irradiated or disrupting planets (e.g., \citealp{brogi2012}, \citealp{rappaport2014}, \citealp{sanchis2015}).

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 \citep{eylen2013}, 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.

\section{Conclusions}

\label{summary}

We have presented a new technique to identify 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 with respect to standard, symmetric models.

We have applied this reasoning to search for hints for clouds in the transits of Kepler-7b, HAT-P-7b and HD209458b.  No clear sign of clouds have been found, although for HAT-P-7b, results seem to indicate a small, yet not statistically very significant, asymmetry consistent with previous phase curve analyses by \citet{vparis2015phase}. This is due to several reasons. 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 HD209458b, 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.

\begin{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 parts by the computing facilities MCIA (M\'esocentre de Calcul Intensif Aquitain) of the Universit\'e de Bordeaux and of the Universit\'e 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.

\end{acknowledgements}

\bibliographystyle{aa}


\section{Introduction}

\mypara{Mixed-criticality Scheduling} In a mixed-criticality system, tasks with different criticality levels
coexist and need to share common resources, such as bandwidth in a communication channel \cite{hts16,hp98}
or execution time on a machine \cite{ves07,bbdlmms10,bls10,cnqw13,pk11,spbb13}.
Contrary to single-criticality systems, the estimated worst-case executing time (WCET) of a task depends on
its criticality level (the higher the criticality level, the more time is estimated). Often, however, the actual
execution times of tasks are not known beforehand and the estimated WCETs deviate hugely from the actual
execution times. The goal is therefore to design {\em robust} schedules that are able to tolerate runtime
variations to a reasonable extent. More conservative WCET estimates are usually used for highly critical tasks
(e.g. braking in a car), while less conservative estimates suffice for low-critical tasks (e.g. displaying the
temperature in a car). In case the allocated time for a task is insufficient at runtime, i.e.,
the actual runtime of a task exceeds its estimated WCET, the execution of the task may nevertheless continue
and suppress subsequent tasks of lower criticality.
For a recent thorough survey on mixed-criticality systems and arising scheduling problems we refer the reader to \cite{bur14}.

\mypara{The Triangle Scheduling Problem}
In this paper, we consider the problem of non-preemptively scheduling $n$ unit length jobs/tasks with different
criticality levels on a single machine. Let $p_i \in \mathbb{N}$ denote the criticality level
of job $i$. The expected execution time of every job is $1$. If however at runtime, a job $i$ requires more time,
then we continue its execution for at most $p_i$ time units, and other jobs with lower criticality levels that were
scheduled in these slots are canceled. The computed schedule has to fulfill the property that, when
prolonging the execution of a job, no other job with larger criticality needs to be canceled. The objective
is to minimize the {\em makespan} (total completion time).

The previous problem can be see as a one-dimensional triangle alignment (or scheduling) problem, defined as follows:

\begin{definition}[Triangle scheduling, Gap] \label{def:ts}
Given integers $p_1, \ldots, p_n$ with $p_1\geq \ldots \geq p_n>0$ find starting times $s_1, \ldots, s_n \geqslant 0$
minimizing the so-called \emph{makespan} $\max_j s_j + p_j$ such that for all $i \neq j$ we have
$| s_i - s_j | \geqslant \min \{ p_i, p_j \}$.
We call \emph{gap} any interval spanned by successive starting times, including the interval between
the maximum starting time and the makespan of the schedule.
\end{definition}
We abbreviate the triangle scheduling problem by \textsc{TS}.

\begin{figure}[ht]

\begin{center}
\begin{tikzpicture}[scale=0.35,thick]
\jobdiscrete{0}{6}{0}
\jobdiscrete{4}{4}{0}
\jobdiscrete{7}{3}{0}
\jobdiscrete{10}{5}{0}
\draw[->](0,0) --  (15,0) node[below] {time};

\draw (0,-2) rectangle (1,-3);
\draw (4,-2) rectangle (5,-3);
\draw (7,-2) rectangle (8,-3);
\draw (10,-2) rectangle (11,-3);
\draw[->](0,-3) --  (15,-3);

\draw (0,-4) rectangle (5,-5);
\draw (7,-4) rectangle (9,-5);
\draw (10,-4) rectangle (14,-5);
\draw[->](0,-5) --  (15,-5);

\draw (0,-6) rectangle (4,-7);
\draw (4,-6) rectangle (8,-7);
\draw (10,-6) rectangle (12,-7);
\draw[->](0,-7) --  (15,-7);

\begin{scope}[shift={(20,0)}]
\job{0}{6}{0}
\job{4}{4}{0}
\job{7}{3}{0}
\job{10}{5}{0}
\draw[->](0,0) --  (15,0) node[below] {time};
\end{scope}
\end{tikzpicture}
\end{center}
  \caption{Example of an optimal schedule. \textbf{Left:} discrete version. \textbf{Below:} different possible executions.  In the first one each job requires a single time slot and all jobs are executed.  In the second and third executions, jobs require longer processing times preventing the execution of some jobs with lower criticality.
  \textbf{Right:} continuous version. Each triangle
  corresponds to a job $j$ labeled with its criticality level $p_j$. The left border of
  a triangle is the starting time of the job, while the right end is its worst-case completion time.}
    \label{fig:ExampleSolution}
\end{figure}

\mypara{Our Results}
In this paper, we initiate the study of \textsc{TS}. We show that the problem is strongly NP-hard,
which implies that \textsc{TS} does not admit a fully polynomial time approximation scheme (FPTAS),
unless $P=NP$. We provide a quasi polynomial time approximation scheme (QPTAS),
which implies that \textsc{TS} is not APX-hard.
In addition we present a Greedy algorithm, which processes the triangles from largest to smallest,
placing each into the largest gap and potentially shifting subsequent triangles to the right
if the gap is not large enough to contain it. We show that this algorithm has an approximation factor
between $1.05$ and $1.5$.

\mypara{Binary Tree Ratio} Furthermore, we establish a measure, denoted the {\em binary tree ratio}, that allows us to
distinguish hard from easy instances:
\begin{definition}[Binary tree ratio]
Given an instance of \textsc{TS} $p=p_1,\ldots,p_n$, we define its \emph{binary tree ratio} $R(p)$ as
\[
R(p) := \max_{i=2, \hdots,n} p_{\lceil i/2 \rceil} / p_i.
\]
\end{definition}
Schedules computed by our Greedy algorithm can be represented by binary trees on the jobs of the
problem instances (see Section~\ref{sec:greedy-tree} details). The binary tree ratio is the maximum ratio between a
job and its successor in the tree. We will show that our Greedy algorithm solves an instance to optimality if its binary
tree ratio is at most $2$. On the other hand, we prove that there are instances
with binary tree ratio strictly larger but arbitrarily close to $2$ that render the problem NP-hard.
A binary tree ratio of $2$ is hence the cut-off point that separates hard from easy instances.

\mypara{Other Related Works}
Since a few decades the scheduling community has been very interested in producing robust schedules that can react to changes in job characteristics.  For example in \cite{HonkompMockus:97:Robust-scheduling} a model is studied where the processing times can vary, and a schedule has to be produced with good objective value even under these variations. For more information see the survey \cite{BriandLa:07:A-robust-approach} as well as a recent PhD thesis and references therein \cite{Wilson:16:Robust-scheduling}.

Vestal \cite{ves07} introduced the \emph{mixed-criticality framework}, where the execution of lower critical jobs can be canceled in order to grant high criticality jobs the necessary amount of resources.  Applications are mostly embedded systems.  In a communication system, jobs represent messages, and safety-critical messages have to co-exist with less critical ones that are not subject to hard constraints.  For instance, the IEC 61508 standard defines four Safety Integrity Levels (SIL) (e.g. the importance of a safety-related job performing braking in a car is much higher than the importance of a job displaying the engine's temperature).
Similarly the \emph{CANaerospace} protocol specifies several criticality levels for messages, and in order to guarantee deliver times of high critical messages, the transmission of lower critical messages can be canceled \cite{Stock:14:CANaerospace-specification}.  This feature has been studied also for the \emph{FlexRay} protocol used in modern cars \cite{dvo14}.

Our model can be seen as a special case of the message transmission model described
by Hanz\'alek et al. in \cite{hts16}. They consider a single machine scheduling problem with release
times, deadlines, different criticality levels, and WCETs that depend on the task and the criticality
level. They propose a linear programming formulation of the problem and prove NP-completeness of this
more general problem. Note that our model does not consider release times and deadlines, and the WCET
of a task equals its criticality level.

Finally we would like to mention a connection with the computational problem of packing triangles in a given rectangle, which has applications in industrial cutting and storage. The later has been shown to be NP-hard~\cite{Chou:16:NP-Hard-Triangle}, while our paper shows that the problem is already hard in the particular case of right triangles that can only be vertically translated and not rotated.

\mypara{Outline}
In Section~\ref{sec:greedy} we consider our Greedy algorithms. Then,
in Section~\ref{sec:NPhardness}, we show that solving the problem for instances with binary tree ratio
strictly larger than $2$ is strongly NP-hard. Last, we provide a quasi polynomial
time approximation scheme (QPTAS) in Section~\ref{sec:QPTAS}.

\section{The Greedy algorithm}
\label{sec:greedy}

We propose a polynomial time approximation algorithm denoted \emph{Greedy}. Recall that
jobs are sorted with respect to their criticality levels, i.e., $p_1 \ge p_2 \ge \dots \ge p_n$.

\begin{definition}[Greedy]
Job 1 starts at time $s_1=0$.  Then, every job $j=2,\ldots,n$, is placed in a largest gap (the first one in case of tie).
 If the chosen gap has length $x$ and starts at time $s_i$, then the current job $j$ is placed at $s_j=s_i+p_j$.  If $2p_j>x$, then all jobs $k$ with $s_k>s_j$ are delayed by $2p_j-x$ in order to maintain feasibility, see Figure~\ref{fig:insert}.
\end{definition}

Note that in case $2p_j>x$ the makespan increases by $2p_j-x$. Hence by choosing the largest gap, Greedy minimizes the increase of the makespan at every step.

\begin{figure}[ht]
\begin{center}
\begin{tikzpicture}[scale=0.18,thick]
\joblabel{0}{10}{}
\joblabel{10}{10}{}
\joblabel{25}{10}{}
\joblabel{35}{10}{}
\joblabel{50}{10}{}
\joblabel{64}{10}{}
\joblabel{29}{4}{$p_j$}
\joblabel{57}{7}{$p_j$}
\draw[<->](0,-1) --  node[auto,swap] {$x$} (10,-1);
\draw[<->](25,-1) --  node[auto,swap] {$p_j$} (29,-1);
\draw[<->](29,-1) --  node[auto,swap] {$x-p_j$} (35,-1);
\draw[<->](50,-1) --  node[auto,swap] {$p_j$} (57,-1);
\draw[<->](57,-1) --  node[auto,swap] {$p_j$} (64,-1);
\end{tikzpicture}
\end{center}
  \caption{When Greedy inserts a job $j$ in a gap of size $x$ (left figure) it creates two gaps. One of size $p_j$ and another either of size $x-p_j$ if $x\geq 2p_j$ (center figure) or of size $p_j$ if $x<2p_j$ (right figure).}
    \label{fig:insert}
\end{figure}

\subsection{Lower bound on the optimum}

A simple lower bound can be obtained by relating gaps to jobs in the schedule
and using the fact that a gap between jobs $i$ and $j$ has size at least $\min\{p_i,p_j\}$.

\begin{lemma}  \label{lem:inner-product-lower-bound}
  Let $S=p_{\lceil n/2\rceil+1}+\ldots+p_n$ be the total processing time of the smaller half of
  the jobs, and $m=p_{(n+1)/2}$ if $n$ is odd and $m=0$ if $n$ is even. The optimal makespan OPT is at least
  \[

      m + 2S.

  \]
\end{lemma}
\begin{proof}
Consider an arbitrary feasible schedule, and let $T$ be its makespan.
To obtain a bound on $T$, we charge every gap to a job as follows:

Map every gap between jobs $i,j$ to the smaller job among them, breaking ties arbitrarily. Map the last gap
between job $j$ and the makespan to job $j$.
Now every job $j$ is the image of $0$, $1$ or $2$ gaps in this mapping. Let $a_j \in \{0, 1, 2\}$ be this number.
There are exactly $n$ gaps, hence we have $\sum_{j=1}^n a_j=n$.
Moreover, since every gap is mapped to a job of no larger size and the total gap size is $T$, we have
\[
     \sum_{j=1}^n a_j p_j \leq T.
\]
The proof follows from the fact that the left hand side is minimized when $a_j=0$ for the larger half of the
jobs, $a_j=2$ for the smaller half of the jobs and possibly $a_{(n+1)/2} = 1$ if $n$ is odd.
\end{proof}
Note that this lower bound can be very weak. Consider a $2$-job instance
with $p_1=M$ and $p_2=1$, for a large $M\geq 2$.  Then the lower bound states that the makespan is at least 2, while
the optimal makespan is $M$.

\clearpage
\subsection{Approximation ratio of Greedy}

\begin{lemma}
  The approximation ratio of Greedy is at most $1.5$.
\end{lemma}
\begin{proof}

Consider the final schedule produced by Greedy on instance $p$. We perform the following transformations:

First, we delay each job by as much as possible, while maintaining feasibility, the makespan and the job order.
This transformation might change the order of the gap sizes, but does not modify the actual gap sizes when viewed as a multi-set.

Second, we define a \emph{truncated} instance $p'$ as follows.
For every job $j$ which starts at some $t$ followed by a gap of size $x$ with $x<p_j$, we set $p'_j=x$.  For all other jobs $j$, we set $p'_j=p_j$, see Figure~\ref{fig:truncated}. Makespan as well as feasibility of the schedule are preserved by the truncation.

\begin{figure}
\begin{center}
\begin{tikzpicture}[scale=0.18,thick]
\job{0}{20}{0}
\job{6}{4}{0}
\job{10}{10}{0}
\job{16}{4}{0}
\job{20}{20}{0}
\job{25}{5}{0}
\job{30}{5}{0}
\job{34}{4}{0}
\job{38}{4}{0}
\draw[dotted] (0,6) -- (6,6);
\draw[dotted] (10,6) -- (16,6);
\draw[dotted] (20,5) -- (25,5);
\draw[dotted] (30,4) -- (34,4);
\end{tikzpicture}
\end{center}
\caption{Transformations on a schedule produced by Greedy: Jobs are delayed (right-shifted) and their sizes truncated (dotted lines).}
\label{fig:truncated}
\end{figure}

We assume that Greedy increased the makespan when placing the last job $n$. This assumption is without loss
of generality, since removing jobs that followed the last makespan increase only decreases the makespan of the optimal schedule,
while the makespan produced by the algorithm is preserved.

We claim that for all $i$ we have $p_n \leq p'_i < 2p_n$.  Indeed when job $n$ was placed, all gaps were of size strictly less than $2p_n$ since the insertion of job $n$ increased the makespan.  Furthermore, by induction, it can be shown that after placing job $j$, all gaps are of size at least $p_j$: By the induction hypothesis, after placing job $j-1$, all gaps were of size at least $p_{j-1}$
which is at least $p_j$, by the assumed ordering of the jobs. Then, no matter how job $j$ is placed, it is impossible that a gap of size smaller than $p_j$ is created (see also Figure~\ref{fig:insert}).

From now on, assume that $n$ is even; the proof for the odd case is similar. Let $A$ be the sum of the larger half of the
sizes among $p'_1,\ldots,p'_n$ and $B$ the sum of the smaller half. The truncation process reduces a job to the size of
gap that follows it. Therefore, the makespan of the schedule produced by Greedy on $p'$  (or on $p$)  is $A+B$.

The previous claim implies that all sizes among $p'_1,\ldots,p'_n$ are within a ratio of two, which implies $A\leq 2B$.

From Lemma~\ref{lem:inner-product-lower-bound} we have $\textrm{OPT}(p')\geq 2B$.
Furthermore, we clearly have $\textrm{OPT}(p')\leq \textrm{OPT}(p)$. We can hence upper bound the makespan of the
schedule produced by Greedy on $p$ as
\[
  A + B \leq 3B \leq \frac32 \textrm{OPT}(p') \leq \frac32 \textrm{OPT}(p).
\]
\end{proof}

Note that this analysis did not use the fact that Greedy places jobs in the largest gap. The crucial property required
in the analysis is the fact that when the placement of a job $j$ increases the makespan, then all gaps are of size strictly
less than $2p_j$.

We were not able to determine the exact approximation factor of Greedy. In Figure~\ref{fig:lowerboundgreedy},
an instance is illustrated that shows that the approximation factor of Greedy is at least $1.05$:
On the instance with jobs of processing times $20,20,10,5,5,4,4,4,4$, Greedy
produces a schedule of makespan $42$, by placing them in the order $20,4,10,4,20,5,5,4,4$, while the optimal schedule
places the jobs in order $20,5,10,5,20,4,4,4,4$ and has makespan $40$.  This example gives the following lower bound.

\begin{lemma}
  The approximation ratio of Greedy is at least $1.05$.
\end{lemma}

We conducted a systematic search for stronger
lower bound constructions, but could only obtain tiny improvements. For example, we
found an instance consisting of $52$ jobs showing a lower bound of $101/96 > 1.052$.

\begin{figure}[ht]
\begin{tikzpicture}[scale=0.16,thick]
\job{0}{20}{0}
\job{5}{5}{0}
\job{10}{10}{0}
\job{15}{5}{0}
\job{20}{20}{0}
\job{24}{4}{0}
\job{28}{4}{0}
\job{32}{4}{0}
\job{36}{4}{0}
\end{tikzpicture}
\hfill
\begin{tikzpicture}[scale=0.16,thick]

\job{0}{20}{0}
\job{4}{4}{0}
\job{10}{10}{0}
\job{14}{4}{0}
\job{20}{20}{0}
\job{25}{5}{0}
\job{30}{5}{0}
\job{34}{4}{0}
\job{38}{4}{0}
\end{tikzpicture}
\caption{The optimal schedule (left) and the schedule produced by Greedy (right) for the lower bound instance.}
\label{fig:lowerboundgreedy}
\end{figure}

\subsection{A case where Greedy is optimal} \label{sec:greedy-tree}

\begin{theorem}
  Greedy is optimal for instances with binary tree ratio at most $2$.
\end{theorem}
\begin{proof}
We show by induction on $j$ that after placing job $j$, there are $2$ gaps of size $p_i$, for every
$\lceil j/2\rceil + 1 \leq i < j$, and either a single gap (if $j$ is odd) or $2$ gaps (if $j$ is even) of size $p_{(j+1)/2}$.  This invariant is true after placing job $1$, where there is a single gap of size $p_1$.  When $j$ is even, the job is placed in the single gap of size $p_{j/2}$.  By the assumption on the binary tree ratio we have $2p_j \geq p_{j/2}$, implying that this gap is replaced by 2 gaps of size $p_j$.
When $j$ is odd, the job is placed in one of the 2 gaps of size $p_{(j-1)/2 + 1}$, and for the same reasons as in the even case the gap is replaced by 2 gaps of size $p_j$.  In both cases the invariant is preserved.

The implication of this observation is that by the lower bound of Lemma~\ref{lem:inner-product-lower-bound}, the schedule produced by Greedy is optimal.
\end{proof}

We can relate the jobs in a tree structure as illustrated in Figure~\ref{fig:tree}. Job $1$ is the root of the tree. Then for every job $j=2,\ldots,n$ inserted into a gap assigned to job $i$, job $j$ is a descendant of job $i$. The result is a single root, connected to a binary tree, which is complete except possibly for the last level, which is left padded.  The job labels on this tree are ordered by levels.  The binary tree ratio of the instance is the maximum ratio between a job and its immediate ancestor in the tree, which was the motivation for the name of this ratio.

\begin{figure}
\tikzset{
  treenode/.style = {align=center, inner sep=0pt, text centered,
    font=\sffamily},
  arn/.style = {treenode, circle, draw=black, text width=1.5em, very thick}
}
\begin{center}
\begin{tikzpicture}[scale=0.7,->,>=stealth',level/.style={sibling distance = 25cm/#1/#1,
  level distance = 1.5cm}]
\node [arn]{1}
child{ node [arn] {2}
    child{ node [arn] {3}
            child{ node [arn] {5}
              child{ node [arn] {9}}
              child{ node [arn] {10}}
            }
            child{ node [arn] {6}
              child{ node [arn] {11}}
              child{ node [arn] {12}}
            }
    }
    child{ node [arn] {4}
            child{ node [arn] {7}
              child{ node [arn] {13}}
              child[missing]
            }
            child{ node [arn] {8}
            }
    }
    }
;
\end{tikzpicture}
\end{center}
\caption{Tree structure of the schedule produced by Greedy on any instance of 13 jobs with binary tree ratio at most $2$.}
\label{fig:tree}
\end{figure}

\section{NP-hardness}
\label{sec:NPhardness}

\begin{theorem}
  \textsc{TS} is strongly NP-hard for instances with binary tree ratio strictly larger than $2$.
\end{theorem}
\begin{proof}
We reduce from the strongly NP-hard numerical $3$-dimensional matching problem, see \cite[problem SP16]{GareyJohnson:79:Computers-and-intractability}.  An instance of this problem consists of integers
\[
      a_1,\ldots,a_n, b_1,\ldots,b_n, c_1,\ldots,c_n, D,
\]
with $D\geq 4$, and for all $i$:
\begin{equation}
      D/4 < a_i,b_i,c_i < D/2.
      \label{eq:NPC-promise-indiv}
\end{equation}
Furthermore, we are guaranteed that
\begin{equation}
 \sum_{i=1}^n a_i + b_i + c_i =  nD.
    \label{eq:NPC-promise-sum}
\end{equation}
The goal is to form $n$ disjoint triplets of the form $(i,j,k)$ with $a_i+b_j+c_k=D$.

Fix some arbitrary large constant $M \geq \frac{5D}{4}$.  The instance consists of $5n$ jobs.
\begin{itemize}
\item There are $n$ jobs $E$ of size $8M+5D$.
\item There are $n$ jobs $F$ of size $4M$.
\item For every $i\in\{1,\ldots,n\}$ there is a job $A_i$ of size $2M+2a_i+D$,
\item as well as a job $B_i$ of size $2M+b_i$,
\item and a job $C_i$ of size $M+c_i+D$.
\end{itemize}

We claim that the instance has a solution of makespan $n(8M+5D)$ if and only if the initial 3-Partition instance has a solution.

For the easy direction, given a solution to the 3-Partition instance we construct a schedule consisting of the concatenation for every triplet $(i,j,k)$ of the jobs $E,A_i,F,B_j,C_k$.  Straightforward verification shows that the resulting schedule has the required makespan.

For the hard direction, consider a solution to \textsc{TS} of makespan $n(8M+5D)$.  Its makespan cannot be smaller, by the presence of the $E$ jobs, that need to be scheduled every $8M+5D$ time units.  They structure the time line into $n$ blocks of equal size $8M+5D$ each.

Suppose that a block contains $k$ jobs of size $x_1\geq \ldots \geq x_k$ plus a single $E$ job.
These jobs create $k+1$ gaps in this block.  Each gap can be assigned to the smaller one among the neighboring jobs, and for an assignment we denote by $a_i\in\{0,1,2\}$ the number of assignments the job of size $x_i$ obtained.

Hence a lower bound on the block size is given by the expression $a_1x_1+\ldots+a_kx_k$ for some $\{0,1,2\}$-weights with $a_1+\ldots+a_k=k+1$.  A valid lower bound is given by $a_1 x_1 + \ldots +a_k x_k$ setting $a_i$ to $2$ for the last $\lceil k/2 \rceil$ indices (corresponding to the smaller jobs), setting $a_{k/2-1}=1$ if $k$ is even and setting $a_i=0$ to the remaining indices.  We will use this lower bound to determine the number of jobs of each type in block.

No block can host two $F$ jobs or more, since $3(4M) > 8M+5D$.  Hence every block contains exactly one $F$ job.

No block can host an $F$ and two $A$ jobs, since (placing weights $a_1=a_2=2$ for the $A$ jobs and $a_3=0$ for the $F$ job)
\[
     4 \left(2M+2\frac D 4  + D\right) = 8M + 6D.
\]
Hence every block contains exactly one $F$ and one $A$ job.

Similarly a block cannot not host an $A$ and $F$ job, together with two $B$ jobs, as setting total $a$-weight $4$ for the $B$ jobs and $1$ for the $A$ job results in the lower bound
\[
   4 \left( 2M+ \frac D 4\right) + 2M + 2\frac D 4  + D =  10M + \frac52 D.
\]
Hence every block contains exactly one $F$ one $A$ and one $B$ job.

Finally a block cannot host an $A$, $F$ and $B$ job, together with two $C$ jobs, as setting total $a$-weight $4$ for the $C$ jobs and $2$ for the $B$ job results in the lower bound
\[
   4 \left(M + \frac D 4 + D\right) + 2 \left( 2 M + 2 \frac D 4 \right)  =  8M + 6D.
\]

In conclusion every block contains an $F$ job, and $A_i$ job, a $B_j$ job and a $C_k$ job with the following lower bound for the space occupied for these jobs, using $a$-weight 2 for jobs $B_j,C_k$ and $a$-weight $1$ for job $A_i$
\[
    2M + 2a_i + D + 2(2M+b_j) + 2(M+c_k + D) = 8M + 3D + 2(a_i+b_j +c_k).
\]
Since this value cannot exceed the size of a block, namely $8M+5D$, every block corresponds to a triplet $(i,j,k)$ with $a_i+b_j+c_k\leq D$.  By the assumption (\ref{eq:NPC-promise-sum}) we have equality for every triplet, and this shows that there is a solution to the 3-Partition instance, see Figure~\ref{fig:block} for illustration.

By construction when the jobs are ordered in decreasing size, they are grouped by types, in the order $F,E,A,B,C$.  Hence the binary tree ratio is determined by the ratio between an $F$ and an $E$ job, between a $B$ and an $F$ job, and between a $C$ and an $A$ jobs.  All these fractions can be made arbitrarily close to $2$ by choosing a large enough value for $M$.
\end{proof}

\begin{figure}[ht]
 \begin{center}%
 \begin{tikzpicture}[scale=0.6,thick]
 \joblabel{0}{8.5}{E}
 \joblabel{2.2}{2.2}{A}
 \joblabel{3.2}{1}{C}
 \joblabel{4.4}{4}{F}
 \joblabel{6.4}{2}{B}
 \end{tikzpicture}
  \caption{A schematic view of a block in the optimal schedule obtained from the reduction.}%
  \label{fig:block}
 \end{center}
\end{figure}

\section{A QPTAS}
\label{sec:QPTAS}

\begin{theorem}
  \textsc{TS} admits a quasi polynomial time approximation scheme.
\end{theorem}
\begin{proof}
The proof starts with a sequence of claims.

\begin{claim}
  Rounding all sizes up to the nearest power of $1 + \epsilon$ changes the optimal makespan
by at most a factor $1 + \epsilon$.
\end{claim}
\begin{proof}
  Take an optimal solution and multiply all processing times and start times by factor $1 + \epsilon$. The solution is still feasible, only the unit has changed. Now round all triangles down to the nearest power of $1 + \epsilon$  while keeping the start times fixed.
\end{proof}

\begin{claim}
We may assume that the ratio $p_1/p_n$ is at most $n/\epsilon$.
\end{claim}
\begin{proof}
  The optimal makespan $\textrm{OPT}$ is at least $p_1$. All triangle with size less than $\epsilon p_1/n$ can be put at the end. We do not need to optimize over these, as they increase the makespan by at most an $\epsilon$ factor.
\end{proof}

From now let the smallest triangle have size $p_n=1$ and the largest have size $p_1 \leq n/\epsilon$. Then $p_1 \leq \textrm{OPT} \leq n p_1 \leq n^2/\epsilon$.

\begin{claim}
We may assume there are only $\lceil \log (n/\epsilon) \rceil +1$  number of different sizes.
\end{claim}
The proof follows from the previous two claims.

\begin{claim}
Restricting start times to values from the set $P=\{0, K, 2K, ..., \lceil n^2/\epsilon \rceil K\}$ for $K = \epsilon p_1/n$ increases the optimal makespan at most by factor $1+\epsilon$.
\end{claim}
\begin{proof}
We modify the start times of the jobs, processing them from left to right. For every job, we move it to the next time in $P$, and move simultaneously all subsequent jobs by the same amount to preserve feasibility.  The value of the solution increases by at most $nK \leq \epsilon p_1 \leq \epsilon \textrm{OPT}$.
\end{proof}
Let $Z$ be the set of all different job sizes after the above rounding.
We have
\begin{align*}
  |Z| & = \lceil \log (n/\epsilon) \rceil +1
  \\
  |P| &\leq \lceil n^2 / \epsilon \rceil
\end{align*}
We design a dynamic programming scheme as follows.
Partition the set of jobs into sets $S$ and its complement $\overline S$. We want to compute the optimal schedule placing first jobs from $\overline S$ and then jobs from $S$.  Only a few parameters from the first schedule influence the possibilities of the second schedule, namely the position of the rightmost triangle of each size. Hence we describe every possible configuration by a vector from $P^Z$, which leads to a quasi polynomial number of configurations. The number of possible subsets $S$ is (in terms of size multiplicities) $O(n^Z)$ which is also quasi-polynomial.

Define $F (C, S)$ to be the optimal makespan of a schedule, placing first $\overline S$ with a configuration $C$ and placing then jobs from $S$.
The goal is to compute $F(e, \{1,\ldots,n\})$, where $e$ is the empty configuration, which for technical reasons assigns to each size $x\in Z$ the starting time $-x$.

The basis cases consist of $F(C,\emptyset)$, where $C$ ranges over all configurations from $P^Z$, some of them might be infeasible.
Then $F (C, S)$
 can be computed from $F (C+j, S-j)$ where $S-j$ is $S$ minus one triangle $j\in S$ and $C+j$ is obtained from $C$ by adding the triangle $j$ to the right of $C$ and placing it as early as possible, i.e.\ at time $\max_{x\in Z} C_x + x$.
 The number of choices for $j$ (in terms of size) is at most $|Z|$. Hence, we need at most $|Z|$ look-ups to compute one value.
\end{proof}

\clearpage
\section{Final remarks}

\shapepar{{3}{0}b{0}\\{0}t{0}{6}\\{6}e{0}}
We introduced a new scheduling problem, motivated by mixed criticality. The novelty lies in its combinatorial structure which defines the contribution to the makespan of each job in a non-local manner. We showed that the problem is strongly NP-hard, thus ruling out the existence of a fully polynomial time approximation schedule.  In addition we provided a quasi polynomial time approximation scheme, ruling out APX-hardness.  Furthermore, we introduced a greedy algorithm for this problem, but still do not understand well its approximation ratio. Closing the gap between the lower bound of $1.05$ and the upper bound of $1.5$ is the main question left open by this paper.
We would like to thank Marek Chrobak and Neil Olver for helpful discussions.
This work is partially supported by PHC VAN GOGH 2015
PROJET 33669TC, the grants FONDECYT 11140566 as well as ANR-15-CE40-0015.