© ESO, 2013

1. Introduction

Observing secondary eclipses of transiting extrasolar planets, where the planet disappears behind its host star, offers a window to study the planetary and orbital properties not attainable by the transit and radial velocity (RV) observations alone. Foremost, eclipses provide an indirect means to study planetary atmospheres. The two main measurables, planet’s albedo and brightness temperature, are strongly coupled with the structure and dynamics of the atmosphere, and can be used to educe information about the physical processes governing the atmosphere. Due to the complexity of the physics involved, this additional empirical knowledge is highly valuable when aiming to understand the planetary atmospheres via theoretical modeling. In addition to the information regarding the atmosphere, eclipses yield information about the planet’s orbit. The eclipse center times and durations allow us to constrain the orbit’s eccentricity and argument of periastron to a higher precision than with RV observations alone. This is especially true for the planets orbiting faint or rapidly rotating stars, for which precise RV observations are difficult to obtain.

A planet’s emergent flux is a combination of reflected stellar light and the planet’s thermal radiation. Both the reflective and thermal properties depend on the structure and composition of the planetary atmosphere. The reflectivity is dominated by the existense or the absence of high-altitude clouds, while the thermal radiation arises mainly from the heating due to the stellar irradiation, from tidal heating (Jackson et al. 2008; Leconte et al. 2010), and from radiogenic heating. Although the reflected light normally dominates at visible wavelengths, the thermal radiation from a hot planet in a short-period orbit may make a major contribution at the red end of the visible spectrum (López-Morales & Seager 2007; Fortney et al. 2008).

The structure of the atmosphere and its reflective properties are tightly coupled, and definitive modeling of the emergent flux is complicated at the best. The composition of the highest optically thick atmosphere layers, and the presence or absence of absorbing gases or reflecting clouds, will determine the reflectivity of a planet.

Broadly, modeling leads us to expect low geometric albedos (A_g < 0.3) at visible wavelengths for the majority of hot Jupiters. Moderately hot (T_eff < 1500 K) cloud-free atmospheres are expected to have small albedos due to the absorption by molecules and alkali metals (Burrows et al. 1997; Marley et al. 1999), while the hottest atmospheres (T_eff > 1500 K) can have high-altitude silicate clouds leading to a high albedo (Sudarsky et al. 2000, 2003). However, strong absorption by TiO and VO gases can dominate over the reflectance from the silicate cloud layer, yielding very low albedos for the hottest atmospheres (Fortney et al. 2008). Such strong absorption has the potential to evoke a stratospheric temperature inversion, which would be observable as anomalously deep eclipses in the infrared. However, the high temperature alone does not assure the existense of a temperature inversion layer. For example, strong stellar UV radiation from active stars can destroy the compounds responsible for the high absorption, thus preventing the birth of an inversion layer (Knutson et al. 2010).

The observed eclipses thus far confirm a complicated picture. Many of the secondary eclipse observations agree with a low albedo (e.g., Charbonneau et al. 1999; Rowe et al. 2008), but recent works by Kipping & Bakos (2011) and Demory et al. (2011) have also shown planets that likely have high albedos.

Recently, Coughlin & López-Morales (2012) carried out a systematic search for eclipses in Kepler Q2 light curves, and note that “majority of the detected planet candidates emit more light than expected owing to thermal blackbody emission in the optical Kepler bandpass, and present a trend of increasing excess emission with decreasing maximum effective planetary temperature”. The authors present possible causes for the trend – including non-local thermodynamic equilibrium or other thermal emission; fluorescent transitions; internal energy generation; and erroneously identifying low-mass stars, brown dwarfs or blends as planets – but do not find a basis to favor any of these solutions over the others.

The CoRoT planets are based on one of the highest-precision data sets available, and all of them have been thoroughly characterized in their respective discovery papers. Here we present the results of a homogeneous search for secondary eclipses in the light curves of all CoRoT planets that have been published until June 2012, namely CoRoT-1b to 23b, including the brown dwarf 15b, but excluding the yet unpublished planet 22b.

CoRoT (Convection, Rotation, and planetary Transits, Baglin et al. 2009) is a space telescope developed jointly by the French space agency CNES in collaboration with ESA, Brazil, and seven European national space agencies. Launched on December 2006, it was the first space telescope dedicated to observing planetary transits. CoRoT offers two time cadences. The survey mode delivers data with a cadence of 512 s, created by stacking 16 exposures of 32 s. For the planet candidates identified during an observing run, a rapid time sampling of 32 s is also available. Furthermore, data of bright (≲R = 13.5 mag) targets have been acquired with three-color photometry, while the fainter ones are observed only in a single passband.

Our analysis is based on Bayesian statistics. The first part of the analysis, eclipse search, uses an approach built on Bayesian model selection (Gregory 2005; Ford et al. 2007). Similar approaches have been successful when searching for planetary signals from RV data (Gregory 2007a,b; Tuomi & Kotiranta 2009; Gregory & Fischer 2010), reanalyzing planetary systems having multiple degenerate solutions (Gregory 2005), and disproving proposed planetary RV signals with insufficient evidence (Tuomi 2011; Gregory 2011). The second part, eclipse characterization, uses Bayesian parameter estimation to derive the posterior distributions for the parameters describing the planet and its orbit (Ford 2005; Gregory 2005).

Two of the CoRoT planets, CoRoT-1b (Barge et al. 2008) and CoRoT-2b (Alonso et al. 2008), have been reported with secondary eclipses in the CoRoT light curves (Snellen et al. 2009; Alonso et al. 2009a,b), that have later been confirmed by Earth-based observations in the infrared (Rogers et al. 2009; Gillon et al. 2009; Alonso et al. 2010; Zhao et al. 2012). These two planets are used to verify our method.

A brief introduction to the Bayesian model selection and parameter estimation approaches follows in Sect. 2. The data used in the analysis are summarized in Sect. 3, and the practical side of our analysis methods in Sect. 4. We present our results in Sect. 5, and discuss the detected eclipse candidates individually in Sect. 6. Finally, we review the main implications of our results and the use of Bayesian model selection in eclipse searches in Sect. 7.

2. Theory

2.1. Bayesian statistics

We start with a short overview on the Bayesian statistics, for an in-depth treatise of the Bayesian approach, see, for example, Gregory (2005), Kass & Raftery (1995), and Robert (2007).

Bayes’ theorem states that the posterior probability distribution P(A|D) for a proposition A with given data D can be derived from the prior probability distribution P(A) and the likelihood of the data for the given proposition, P(D|A), as $\mathit{P}\mathrm{\left(}\mathit{A}\mathrm{\right|}\mathit{D}\mathrm{)}\mathrm{=}\frac{\mathit{P}\mathrm{\left(}\mathit{A}\mathrm{\right)}\mathit{P}\mathrm{\left(}\mathit{D}\mathrm{\right|}\mathit{A}\mathrm{)}}{\mathrm{\int }\mathit{P}\mathrm{\left(}\mathit{A}\mathrm{\right)}\mathit{P}\mathrm{\left(}\mathit{D}\mathrm{\right|}\mathit{A}\mathrm{)}\mathrm{d}\mathit{A}}\mathrm{·}$(1)For a model selection problem, different propositions correspond to different models, while for parameter estimation the propositions correspond to different values of a model-specific parameter vector θ.

If we assume the observed light curve fluxes F_O to follow from a model M(θ,t), where t is the exposure center time, and further assume independent and identically distributed (i.i.d.) errors (deviations from the model) that follow a normal distribution with a standard deviation σ, the probability P(e_i|θ) for an error e_i = F_O(t_i) − F_M(θ,t_i) of the light curve point i is $\mathit{P}\mathrm{\left(}{\mathit{e}}_{\mathit{i}}\mathrm{\right|}{\theta }\mathrm{)}\mathrm{=}\frac{\mathrm{1}}{\mathit{\sigma }\sqrt{\mathrm{2}\mathit{\pi }}}\exp \left(\mathrm{-}\frac{{\mathit{e}}_{\mathit{i}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}\right)\mathit{,}$(2)and the likelihood (joint probability) for N i.i.d. errors is obtained using the multiplication rule of independent probabilities $\begin{array}{ccc}\mathit{P}\mathrm{\left(}\mathit{D}\mathrm{\right|}{\theta }\mathrm{)}& & \mathrm{=}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{N}}{\mathrm{􏽙}}}\frac{\mathrm{1}}{\mathit{\sigma }\sqrt{\mathrm{2}\mathit{\pi }}}\exp \left(\mathrm{-}\frac{{\mathit{e}}_{\mathit{i}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}\right)\mathit{,}\\ & & \mathrm{=}{\left(\mathit{\sigma }\sqrt{\mathrm{2}\mathit{\pi }}\right)}^{\mathrm{-}\mathit{N}}\exp \left(\mathrm{-}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}\frac{{\mathit{e}}_{\mathit{i}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}\right)\mathit{,}\\ & & \mathrm{=}{\left(\mathit{\sigma }\sqrt{\mathrm{2}\mathit{\pi }}\right)}^{\mathrm{-}\mathit{N}}\exp \mathrm{(}\mathrm{-}{\mathit{\chi }}^{\mathrm{2}}\mathit{/}{\mathrm{2}}^{\mathrm{)}}\mathit{.}\end{array}$(3)Since the number of points in the CoRoT light curves is large, it is necessary to carry out the calculations using the natural logarithm of the likelihood, which is now $\ln \mathit{P}\mathrm{\left(}\mathit{D}\mathrm{\right|}{\theta }\mathrm{)}\mathrm{=}\mathrm{-}\frac{\mathit{N}}{\mathrm{2}}\ln \mathrm{2}\mathit{\pi }\mathrm{-}\mathit{N}\ln \mathit{\sigma }\mathrm{-}{\mathit{\chi }}^{\mathrm{2}}\mathit{.}$(4)Finally, the logarithm of the posterior probability for a model M with a parameter vector θ given the observed flux F_O is $\ln {\mathit{P}}_{\mathit{M}}\mathrm{\left(}{\theta }\mathrm{\right|}{\mathit{F}}_{\mathrm{O}}\mathrm{)}\mathrm{=}\ln \mathit{P}\mathrm{(}{\theta }\mathrm{)}\mathrm{-}\frac{\mathit{N}}{\mathrm{2}}\ln \mathrm{2}\mathit{\pi }\mathrm{-}\mathit{N}\ln \mathit{\sigma }\mathrm{-}{\mathit{\chi }}^{\mathrm{2}}\mathit{.}$(5)

2.2. Parameter estimation

Parameter estimation in its simplest form boils down to finding the posterior maximum for the model-specific parameter vector. If the posterior is likelihood dominated, this equals to maximizing the likelihood function, and, if the error is treated as constant, minimizing χ².

For a more robust analysis, we usually want to compute the full posterior distributions for the parameters, from which we can derive different summary statistics. In the current study we compute the final posterior distributions using Markov chain Monte Carlo (MCMC) sampling.

2.3. Model selection

In order to compare the ability of our competing models to explain the data, we need to calculate the global likelihood for a model by marginalizing (integrating) the posterior probability over the model parameters as $\mathit{P}\mathrm{\left(}\mathit{M}\mathrm{\right)}\mathrm{=}\mathrm{\int }\mathit{P}\mathrm{\left(}{\theta }\mathrm{\right)}\mathit{P}\mathrm{\left(}\mathit{D}\mathrm{\right|}{\theta }\mathrm{)}\mathrm{d}{\theta }\mathit{.}$(6)Given the global likelihoods computed for each model, we can calculate the Bayes factor (Kass & Raftery 1995; Gregory 2005; Jeffreys 1961) of a model M_i over model M_j as ${\mathit{B}}_{\mathit{ij}}\mathrm{=}\mathit{P}\mathrm{\left(}{\mathit{M}}_{\mathit{i}}\mathrm{\right)}\mathit{/}\mathit{P}\mathrm{\left(}{\mathit{M}}_{\mathit{j}}\mathrm{\right)}\mathit{,}$(7)and given the Bayes factors, we can calculate the normalized model probabilities for all the models as $\mathit{P}\mathrm{\left(}{\mathit{M}}_{\mathit{i}}\mathrm{\right)}\mathrm{=}\frac{{\mathit{B}}_{\mathit{i}\mathrm{0}}}{{\sum }_{\mathit{j}\mathrm{=}\mathrm{0}}^{\mathit{N}}{\mathit{B}}_{\mathit{j}\mathrm{0}}}\mathrm{·}$(8)With two models, the probability for the model M₁ is given by $\mathit{P}\mathrm{\left(}{\mathit{M}}_{\mathrm{1}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}\mathrm{1}\mathit{/}{\mathit{B}}_{\mathrm{10}}}\mathit{,}$(9)and for the model M₀ we get P(M₀) = 1 − P(M₁).

The Bayes factors can also used directly to give “a rough descriptive statement about standards of evidence” (Kass & Raftery 1995), and we reproduce the guidelines for the direct interperation as given by Kass & Raftery in Table 1.

2.4. The role of prior distributions

The prior distributions (priors) encode our knowledge about the plausible parameter values and relative model probabilities. In parameter estimation, the importance of a prior depends on the relative amounts of parameter-specific information encoded by the prior and information incorporated into the data. If we already have an estimate for a parameter, e.g., from previous research, we can use the estimate to construct an informative prior that may have a substantial role in the posterior distribution. In the absence of a prior estimate, we use an uninformative prior that aims to maximize our uncertainty regarding the parameter value, and the posterior will be dominated by the likelihood.

We use two uninformative priors in our work. The Uniform prior, $\mathit{P}\mathrm{\left(}\mathit{x}\mathrm{;}\mathit{a,b}\mathrm{\right)}\mathrm{=}\mathrm{1}\mathit{/}\mathrm{(}\mathit{b}\mathrm{-}\mathit{a}\mathrm{)}\mathit{,} \mathit{a}\mathrm{\le }\mathit{x}\mathrm{\le }\mathit{b,}$(10)is used for location parameters, and the Jeffreys’ prior, $\mathit{P}\mathrm{\left(}\mathit{x}\mathrm{;}\mathit{a,b}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathit{x}\ln \mathrm{(}\mathit{b}\mathit{/}\mathit{a}\mathrm{)}}\mathit{,} \mathrm{0}\mathit{<}\mathit{a}\mathrm{\le }\mathit{x}\mathrm{\le }\mathit{b,}$(11)is used for scale parameters (for a more detailed explanation of location and scale parameters, see, e.g., Gregory 2005). The parameters a and b define the span of the prior. For the MCMC characterization we also use a Normal prior $\mathit{P}\mathrm{\left(}\mathit{x}\mathrm{;}\mathit{\mu ,\sigma }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathit{\sigma }\sqrt{\mathrm{2}\mathit{\pi }}}\exp \left(\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{\mathit{x}\mathrm{-}\mathit{\mu }}{\mathit{\sigma }}\right)}^{\mathrm{2}}\right)\mathit{,}$(12)where μ is the mean and σ is the standard deviation of the distribution.

As long as the prior distribution encompasses the likelihood function, the exact form of an uninformative prior has often only a small effect in parameter estimation. However, in model selection, the prior widths reflecting the state of our ignorance can affect the global likelihood significantly.

The greater our uncertainty about a parameter value is, the wider the prior will be. Since priors are normalized probability distributions, a wider prior yields a smaller posterior point value, and a smaller marginal likelihood.

Thus, if our prior distributions are wider than the likelihood distributions, each additional model parameter adds a penalization factor to the global likelihood. These factors are multiplicative, and work to penalize complex models against simpler ones. For a complex model to be considered better than a simpler one, the additional complexity must be supported by a sufficiently improved likelihood.

3. Data

For the analysis of CoRoT-1b to 21b and 23b, we use the latest versions of the CoRoT N2 light curves available publicly from the IAS Data Center¹ as of June 2012. We use the chromatic light curves when available, but include from them only the red channel for the analysis to reduce the number of jumps from cosmic ray hits. We use both the 32 and 512 second time cadences (when available) by assigning the two cadences a separate mean point-to-point scatter (error) estimate in the analysis, and exclude the phase-span near the primary transit completely from the analysis.

The stellar and planetary properties used in constructing the priors are gathered from the latest publications discussing the systems. The adopted planet properties are listed in Table 2 and the host star properties in Table 3, both sorted by the presence of eclipse events, as indicated in Sect. 6.

4. Analysis

4.1. Overview

Our analysis consists of three main steps. First, we identify the eclipse candidates and assess their significances. Next, we carry out tests for the plausible candidates to verify that they are not produced by singular events (jumps) or correlated noise. Finally, we characterize the plausible eclipse candidates that pass our tests by calculating the marginal posterior probabilities for the eclipse model parameters.

Our code combines Fortran2003 for the numerically intensive computations and Python² for the high-level functionality. We use simple Monte Carlo (MC) importance sampling to obtain the global likelihood estimates, and MCMC methods to obtain the parameter posterior distributions. Parallelization is carried out using OpenMP and MPI. We use the NumPy³, SciPy³, MPI4Py⁴, matplotlib⁵, PyFITS⁶, and emcee (Foreman-Mackey et al. 2012) – libraries for the data analysis, MPI parallelization, and data IO.

We use the Mandel & Agol (2002) transit shape model without limb darkening to model the eclipse. The model is implemented in Fortran2003 and uses OpenMP for parallelization within a single computing node.

The computations were carried out using the Diodo cluster located at the Instituto de Astrofísica de Canarias and the LaPalma supercomputer located at the La Palma Astrophysics Centre.

4.2. Models and parameterization

We consider two models: with and without an eclipse (M₁ and M₀, respectively). Table 4 offers an overview of the model parameters and their priors. With M₀ being a submodel of M₁, the models share parameters, which allows us to sample both of the model likelihood spaces simultaneously.

The no-eclipse model, M₀, attempts to explain the data with a constant zeropoint value, θ_f0, and mean point-to-point scatter estimates for the light curves of both time cadences (θ_σl and θ_σs, respectively). Thus, M₀ has 2–3 free parameters, depending on whether both time samplings are available for a given CoRoT target.

The eclipse model, M₁, includes the parameters of M₀, and introduces eight additional parameters. However, only three of these additional parameters have posterior distributions dominated by the likelihood from the data. The five other additional parameters have prior-dominated posteriors, and are included to propagate the uncertainties in their estimates into the analysis.

Of the prior-dominated parameters, the radius ratio, period, scaled semi-major axis, and inclination have priors based on published values that were estimated from the primary transits. The contamination factor gives the amount of third light inside the CoRoT aperture mask (Deeg et al. 2009; Deleuil et al. 2009), and has a prior based on the published contamination values.

Of the likelihood-dominated parameters, the planet-star flux ratio θ_f gives the ratio of the surface brightness of the planet to the surface brightness of the star, and is related to the eclipse depth ΔF and planet-to-star surface area ratio ${\mathit{\theta }}_{\mathrm{k}}^{\mathrm{2}}$ as ${\mathit{\theta }}_{\mathrm{f}}\mathrm{=}\mathrm{\Delta }\mathit{F}\mathit{/}{\mathit{\theta }}_{\mathrm{k}}^{\mathrm{2}}$. The eccentricity and argument of periastron determine the eclipse center and duration. Their posterior probabilities are constrained by radial velocity orbservations, but the constraints are often wide.

Both models also include two derived parameters: the eclipse center and duration, T_e and D_e, respectively. The dependency of M₀ on these two parameters is due to our chosen detrending approach, and is explained in more detail in Sect. 4.4.

4.3. Parameter priors

As mentioned previously, priors have an important role in the Bayesian model selection. Since the eclipse model, M₁, has eight parameters more than the no-eclipse model, M₀, it would be expected that the penalization against M₁ would be strong.

However, the five prior-dominated M₁ parameters are closer to fixed constants than to free parameters, and the penalization by most of them has a relatively small effect (the likelihood within the prior limits is close to a constant). For these parameters, the additional information in the eclipse is very small compared to the information from the primary transit and radial velocity observations. The main motivation for including the prior-dominated parameters comes from the added robustness. The parameters are included as nuisance parameters to be marginalized over in order to propagate their uncertainties to the global likelihoods.

Excluding the unknown flux ratio θ_f, all the M₁ parameters have uniform priors centered around the latest published values and half-widths corresponding to their 1σ uncertainties. In many cases the eccentricities are poorly constrained by the radial velocity measurements, and we set the priors to span the eccentricity space from zero to the published maximum value.

The zeropoint and point-to-point scatter are handled somewhat differently than the other parameters. The posteriors of these three parameters were found to be close to normal when the other parameters are held fixed. This allows us to use Laplace’s approximation (see, for example, Gregory 2005) to obtain an analytical likelihood estimate marginalized over the three parameters for each MC sample, and reduce the sampling variance without sacrificing accuracy.

4.4. Light curve detrending

Since the depth of an eclipse is on the order of percents of the depth of a transit, great care must be taken to remove the signals not related to an eclipse. CoRoT light curves contain four major sources of correlated noise: stellar variability on long time-scales, short time-scale variability from stellar granulation (Aigrain et al. 2009), random jumps due to cosmic ray hits (Pinheiro da Silva et al. 2008), and a quasi-periodic signal related to the satellite’s orbit (Pinheiro da Silva et al. 2008; Aigrain et al. 2009).

We chose a detrending approach in which the detrending is carried out separately for each planetary period, and depends on the period, eclipse center and eclipse duration. Thus, we need to detrend the data separately for each MC integration sample. This increases the computational burden (the detrending is one of the most computationally intensive tasks in the analysis), but allows us to average over the possible effects caused by the detrending.

The primary source of light curve variability is usually the star itself. Fortunately, since the time scales for stellar variability are significantly longer than the eclipse duration, this variability is also the easiest external signal to remove: the trends can be modeled locally with a polynomial fitted to the out-of-eclipse (OE) data around each eclipse.

The second major noise signal is due to high-energy particles impacting the detector (Pinheiro da Silva et al. 2008). The impacts cause sudden jumps in the flux level, followed by either gradual or abrupt return close to the original level. While methods have been developed to correct for these jumps (Mislis et al. 2010), we decide to simply exclude the orbital periods with identifiable jumps within the eclipse search range. This is a preferable approach since even a small error in the correction could still yield a signal strong enough to mimic an eclipse.

The third noise signal comes from the crossing of the South Atlantic Anomaly (SAA) and the satellite’s entering and exiting from the Earth’s shadow (Pinheiro da Silva et al. 2008; Aigrain et al. 2009). The signal is quasi-periodic over the period of the satellite’s orbit, has relatively small amplitude, and does not change significantly over the duration of an eclipse. We first reduce the signal by excluding the points flagged by the N2 pipeline as having been obtained during the SAA crossing. The residual signal is modeled and removed using a periodic spline fitted to the OE flux around each eclipse folded over the satellite’s period.

Finally, while the variability due to stellar granulation has a low amplitude, it can still be orders of magnitude greater than an eclipse signal (especially, K-dwarfs and giants have been shown to feature variations on a 0.5 mmag level on 2 h time scales, Aigrain et al. 2009). Since the variability is on the same time-scale as the eclipses, it cannot be removed safely. The Bayes factor mapping of light curves with substantial short time-scale variation usually shows several high-probability peaks with unrealistic eclipse depths. These peaks from single correlated-noise events can often be identified at the later stages of the analysis, as explained in Sect. 4.6, but may also render the analysis of faint targets infeasible.

Fig. 1 Main results for CoRoT-1b. Shown are the eclipse model probability and the Bayes factor (a) and b), respectively); the lightcurve and the model, both folded over the period and binned with three different bin widths c); and posterior distributions for the flux ratio, eclipse center phase, eccentricity and the argument of the periastron (d)–g)).

4.5. Eclipse search

The eclipse search is based on Bayesian model selection. We carry out the most numerically intensive task, the likelihood sampling, separately from the rest of the analysis. After we have the likelihood samples, the computation of posterior probabilities, marginal posteriors, global likelihoods and Bayes factors is fast. Also, this separation allows us to evaluate the effect that different priors have on the posteriors and on the Bayes factors.

We start by calculating a set of likelihood samples for the two competing models, M₀ and M₁, using MC importance sampling combined with the Laplace approximation (LA). The MC sampling ranges are based on the published parameter values with their errors. The Laplace approximation for the likelihood marginalized over the zeropoint and point-to-point scatter is calculated for each MC sample. The approximation is calculated by first maximizing the likelihood as a function of the LA parameters, then computing the covariance matrix in these dimensions around the maximum, and finally using the obtained covariances to calculate the integral (for a detailed example, see Gregory 2005, p. 291).

Given the likelihood samples for both models, we calculate the Bayes factors for a set of eclipse centers, C. The set spans the whole ecplice center space (determined by the θ_p, θ_e, and θ_ω priors) linearly. We compute the global likelihoods with a constant-width uniform prior on the eclipse center time. For each T_e ∈ C, the prior is centered around T_e, and has a width of one tenth of the transit duration. The method generates a one-dimensional Bayes factor map, from which the eclipse candidates can be identified (for example, see Fig. 1). The approach corresponds to making n_C model comparisons, where n_C is the number of suggested eclipse centers, where each comparison differs by the priors imposed on the parameters.

4.6. Plausibility tests

After the initial identification of the eclipse candidates, further steps must be taken to assess their significance. The most important criterion we use for testing the plausibility of a candidate is the “splitting test”, where we require that the signal must be found from separate subsections of the light curve. A jump by a high-energy particle hit or correlated noise, e.g., from stellar activity, can both generate dips that mimic an eclipse signal, but these signals arise from single events.

To separate signals arising from single events from plausible eclipse signals, we map the Bayes factors for each CoRoT target featuring a significant signal for three light curve subsets. We remove one thirds of the total light curve from each subset, and require that the eclipse signal is visible for all the subsets. This approach allows us to reject most of the false eclipse candidates arising from single events, but cannot be used to test the weakest signals.

4.7. Eclipse characterization

After identifying the most plausible eclipse candidates, we carry out MCMC simulations to obtain the posterior distributions for the M₁ parameters. We originally used our own implementation of the Metropolis-Hastings sampler, but changed to use emcee by Foreman-Mackey et al. (2012), a freely available Python implementation of the affine invariant Markov chain sampler with excellent convergence properties (Goodman & Weare 2010). The change of the sampler did not alter the results, but reduced the computation time significantly due to the improved posterior sampling properties (smaller number of MCMC steps were needed in order to obtain independent samples).

We change our parameterization slightly for the characterization step. First, we assign proper priors for the zeropoints and errors (uniform and Jeffreys’, respectively). Second, since we already have an estimate for the eclipse center from the Bayes factor mapping, T_max(B10), we assign a normal prior N(μ = T_max(B10), σ = 0.2 h) on the eclipse center. We also add a prior on the primary transit duration, which allows us to constrain the (θ_e, θ_ω)-space.

5. Results

We summarize our results for the planets with eclipse candidates in Table 5 and Figs. 1 to 5. The first two columns in Table 5 list the Bayes factors and the corresponding M₁ posterior probabilities from the model comparison analysis. The next four columns show the median estimates for the eclipse center, eccentricity, flux ratio and eclipse depth from the MCMC analysis. The uncertainties correspond to the 68.2% credible intervals, or the 1σ confidence limits if the posterior distributions are close to normal. Finally, the last two columns list the 95% credible intervals for the brightness temperature (the temperature estimated from the eclipse depth assuming the planet to radiate as a black body), T_br, and equilibrium temperature (the temperature the planet would have under thermal equilibrium if the stellar radiation would be the only heat source), T_eq.

The panels in the Figs. 1–5 plot (a) the model probability and (b) Bayes factor as a function of the eclipse center time; (c) the CoRoT light curve and the model derived from the MCMC analysis, both phase-folded and binned using three bin widths; and (d) the marginal posterior distributions from the MCMC analysis for the flux ratio, (e) eclipse center phase, (f) orbit eccentricity, and (g) argument of the periastron.

Fig. 2 As in Fig. 1, for CoRoT-2b.

We list the secondary eclipse ephemerides in Table 6, and show the expected flux ratio, eclipse depth and equilibrium temperature limits for planets without eclipse candidates in Table 7. The equilibrium temperatures are calculated assuming the planet to radiate as a black body, and the estimates are calculated using the published values for the effective stellar temperature, planet radius ratio and semi-major axis.

The brightness and equilibrium temperature estimates are calculated in Monte Carlo fashion combining the published estimates for the semi-major axis, radius ratio, and stellar temperature with our flux-ratio estimate. The equilibrium temperatures are calculated as ${\mathit{T}}_{\mathrm{eq}}\mathrm{=}{\mathit{T}}_{\mathit{\star }}{\mathit{\theta }}_{\mathrm{a}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}{\left(\mathit{f}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{A}}_{\mathrm{B}}\mathrm{)}\right)}^{\mathrm{1}\mathit{/}\mathrm{4}}\mathit{,}$(13)where T _⋆  is the estimated stellar temperature, θ_a the scaled semi-major axis, f the heat redistribution factor, and A_B the Bond albedo. The brightness temperature, T_br, is numerically solved from the equation ${\mathit{\theta }}_{\mathrm{f}}\mathrm{=}{\mathit{A}}_{\mathrm{g}}{\mathit{\theta }}_{\mathrm{a}}^{\mathrm{2}}\mathrm{+}\mathit{B}\mathrm{\left(}{\mathit{T}}_{\mathrm{br}}\mathrm{\right)}\mathit{/}\mathit{B}\mathrm{\left(}{\mathit{T}}_{\mathit{\star }}\mathrm{\right)}\mathit{,}$(14)where the first term represents the reflected light and the second the thermal radiation with θ_f being the surface flux ratio, A_g the geometric albedo (we use A_g = 1.5A_B), and B is Planck’s law (i.e., we approximate the stellar radiation with a blackbody). The planet’s Bond albedo is allowed to range from 0 to 0.3, and the heat redistribution factor from 1/4 to 2/3. Thus, the reported 95% credible intervals advises us on the possible values of T_br and T_eq assuming a small albedo and without assumptions about the heat redistribution efficiency.

Fig. 3 As in Fig. 1, for CoRoT-6b.

Fig. 4 As in Fig. 1, for CoRoT-11b.

Fig. 5 As in Fig. 1, for CoRoT-15b.

6. Discussion

6.1. Non detections

The Bayes factor maxima in support of the model M₁ for the planets listed in Table 7 are all below unity, or the peaks do not pass the plausibility tests. In most cases the planets without detected eclipse candidates have a small expected T_eq, their host star shows strong small-scale variability, or the host star is relatively faint.

The only major exception is CoRoT-7b with a high equilibrium temperature and a bright host star. However, the small size of CoRoT-7b makes the detection of its eclipse highly unlikely.

CoRoT-14b has the highest expected equilibrium temperature of the set, but its host star is very faint (V = 16), and the light curve is dominated by correlated small time-scale noise. The Bayes factor mapping of CoRoT-14b light curve shows several peaks, but they can all be attributed to single events due to correlated noise.

6.2. Previously reported eclipses

6.2.1. CoRoT-1b

We use CoRoT-1b (Barge et al. 2008) as a first test case for our method. The light curve has been shown to contain an eclipse (Alonso et al. 2009b; Snellen et al. 2009), which has been confirmed with ground-based near-infrared (NIR) observations by Rogers et al. (2009), Gillon et al. (2009) and Zhao et al. (2012).

We identify the eclipse with an M₁ posterior probability of 88%. The results from the model comparison and parameter estimation are presented in Table 5 and Fig. 1. Our flux ratio estimate of 1% ± 0.4% yields an eclipse depth of (2 ± 0.8) × 10^-4, which is within the 1σ limits of the white light curve depth of (1.6 ± 0.6) × 10^-4 reported by Alonso et al. and the red channel depth of (1.3 ± 0.4) × 10^-4 obtained by Snellen et al.

Compared to the monochromatic depth by Alonso et al., the marginally deeper transit in the red channel is expected if the planet has a small albedo in visible light. Our results diverge more strongly from those by Snellen et al., who find a smaller depth in the red channel than Alonso et al. in the white light. However, the differences are not statistically significant, and can be due to the differences in the data reduction.

6.2.2. CoRoT-2b

CoRoT-2b (Alonso et al. 2008) is used as a second test case. The light curve contains an eclipse signal (Alonso et al. 2009b) that has again later been confirmed by ground-based NIR observations (Alonso et al. 2010; Cacéres et al., based on VLT observations, pers. comm.).

We identify the eclipse with a marginal M₁ probability of 66% (Bayes factor of 2). Alonso et al. (2009b) report a secondary eclipse depth of (6.6 ± 2) × 10^-5 at orbital phase 0.494 ± 0.006, values that agree with our estimates of a depth of (7 ± 5) × 10^-5 at phase 0.501 ± 0.002. Our uncertainty estimate for the depth is considerably greater than the estimate by Alonso et al., which is probably due to the different approaches used for the parameter estimation.

6.3. Statistically significant eclipse candidates

6.3.1. CoRoT-6b

CoRoT-6b is a hot Jupiter with M_P = 2.96 M_Jup and R_P = 1.17 R_Jup orbiting an F9V star on a relatively long-period orbit of 8.9 days (Fridlund et al. 2010). The light curve features an eclipse candidate with a small but non-zero eccentricity (θ_e = 0.06 ± 0.01). We obtain a Bayes factor of 10.6 (see Fig. 3), corresponding to posterior M₁ probability of 91.4%. The candidate passes all our tests, and does not seem to be due to a single event.

However, the candidate is deep (θ_f = 1.96% ± 0.69%) considering the planet’s period, and yields a minimum planetary brightness temperature that is nearly 1000 K higher than the maximum equilibrium temperature. If the eclipse is real, substantial amount of additional heating would be needed to explain the depth.

6.3.2. CoRoT-11b

The light curve of CoRoT-11b (Gandolfi et al. 2010) shows the most significant eclipse candidate identified by our study. We obtain a Bayes factor of (5.5 ± 1.1) × 10⁴, corresponding to a M₁ probability of 99.998%, for an eccentric orbit with θ_e = 0.35 and θ_ω = 75°. The candidate signal passes the tests, and is unlikely to be due to a single event or correlated noise.

CoRoT-11b is an inflated planet with M_P = 2.33 M_Jup and R_P = 1.4 R_Jup, orbiting on a 2.99 d period a hot rapidly rotating F6V star. Gandolfi et al. (2012) measure a sky-projected spin-orbit angle of 0.1° ± 2.6° using Doppler tomography, and conclude that the orbit of the planet is most likely closely aligned with the stellar rotation axis.

The rapid rotation of the primary makes accurate estimation of the eccentricity from radial velocity observations difficult. Gandolfi et al. (2010) give an RV-based constraint of θ_e < 0.7, and comment that for θ_e ≳ 0.2 the mean stellar density ρ estimated from the transit would be incompatible with an F6 dwarf star (the ρ reported by Gandolfi et al. is 0.69 ± 0.02 g/cm³). Using the Eq. (4) for the mean stellar density⁷ρ given by Tingley et al. (2011a), $\mathit{\rho }\mathrm{=}\frac{\mathrm{3}{\mathit{\theta }}_{\mathrm{p}}{\mathit{Q}}^{\mathrm{3}}}{{\mathit{\pi }}^{\mathrm{2}}\mathit{G}{\mathit{\tau }}_{\mathrm{14}}^{\mathrm{3}}}{\left(\frac{\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{\theta }}_{\mathrm{k}}{\mathrm{)}}^{\mathrm{2}}\mathrm{-}{\mathit{b}}^{\mathrm{2}}{\mathit{Q}}^{\mathrm{2}}}{\mathrm{1}\mathrm{-}{\mathit{\theta }}_{\mathrm{e}}^{\mathrm{2}}}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathit{,}$(15)where τ₁₄ is the transit duration, G the gravitational constant, and $\mathit{Q}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\theta }}_{\mathrm{e}}^{\mathrm{2}}\mathrm{)}\mathit{/}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{\theta }}_{\mathrm{e}}\sin {\mathit{\theta }}_{\mathit{\omega }}\mathrm{)}$, with our θ_e and θ_ω values, we obtain ρ = 0.51 g/cm³. Gandolfi et al. state that ρ would be unrealistically high for high eccentricities, but this applies only for θ_ω = 0. However, if we decrease the impact parameter to b = 0.39 (necessary if the eclipse is considered to be real, see below), we can reproduce the reported value of ρ = 0.69 ± 0.02 g/cm³, and ρ is within 3σ from its reported value for 0.15 < b < 0.55.

The value of θ_ω = 75° implies that the eclipse occurs near the apoastron, and the transit near the periastron. The impact parameter for the transit must therefore be smaller then the impact parameter for the eclipse. The current transit-derived impact parameter estimates of 0.8 (Gandolfi et al. 2010, 2012; Southworth 2011) have been calculated assuming a circular orbit, and only a dedicated light curve modeling beyond the scope of this work could show if the transit could be reproduced with the given θ_e, θ_ω and a close-to-zero impact parameter.

We derive a planet-to-star flux ratio of 3.2%, leading to an eclipse depth of (3.6 ± 0.7) × 10^-4. This corresponds to T_br = 2800 K, which is significantly higher than the maximum equilibrium temperature of 2200 K. The high temperature might be explained by tidal heating due to high eccentricity (Jackson et al. 2008; Leconte et al. 2010), but the analysis of the physical plausibilities of different heating mechanisms is again outside the scope of this work. Further RV measurements to better constrain the eccentricity are subject to ongoing observations.

Given the deep eclipse in the visible light and the brightness of the primary (V = 12.8), the eclipse should be confirmable using ground-based NIR observations. Rough predictions for the J, H, and K band depths are 1.0, 1.6, and 2.2 mmags.

6.3.3. CoRoT-15b

CoRoT-15b (Bouchy et al. 2010) is an inflated high-mass brown dwarf with M = 63 M_Jup and R = 1.1 R_Jup on a P = 3 day orbit around an F7V star. The star is faint (V ≈ 16), and the light curve is monochromatic, contains only long-cadence data, and shows several jumps and strong stellar variability.

We identify an eclipse candidate with a Bayes factor of 27, corresponding to a M₁ probability of 97%. The candidate signal is deep, yielding a minimum brightness temperature estimate of 3500 K. The temperature is suspiciously high, but given the high mass, inflated radius, stellar insolation, and possibly young age of the system (Bouchy et al. give two age estimates derived using different methods, 1.14–3.35 Gyr and 1.9 ± 1.7 Gyr, where the latter method gave possible pre-main sequence solutions), may still be physically feasible.

6.4. Marginal detections

In addition to the statistically significant candidates, we also report four marginal candidates that pass the plausibility tests. Considering their Bayes factors, all the marginal candidates fall into the “not worth more than a bare mention” class (Jeffreys 1961; Kass & Raftery 1995).

1. CoRoT-3b is an object near the planet/brown dwarf boundary (M_P = 21.66 M_Jup, R_P = 1.01 R_Jup) in a 4.26 day orbit (Deleuil et al. 2008). We find a tentative eclipse near an orbital phase of 0.51. The candidate is marginal with a Bayes factor B₁₀ of 2.3, and the light curve shows correlated noise on the time-scale of hours. The 95% credible intervals for T_br and T_eq overlap, but the median brightness temperature (2700 K) is high compared to the upper limit equilibrium temperature of 2200 K. However, with its large mass, CoRoT-3b is above the deuterium-burning mass limit (Spiegel et al. 2011), and the additional heating from deuterium burning early during the planet’s history might be able to explain the temperature.

2. CoRoT-13b is a high-density hot Jupiter with M_P = 1.3 M_Jup, R_P = 0.9 R_Jup, and P = 4.04 days (Cabrera et al. 2010). The light curve features several jumps, and the presence of low-amplitude small-time-scale variability is evident. We find a tentative eclipse passing our tests with B₁₀ = 3.3 near a phase of 0.48. Unlike CoRoT-3b, the T_br and T_eq 95% credible intervals do not overlap, and the median T_br (2800 K) is significantly higher than the upper limit equilibrium temperature (1700 K).

3. CoRoT-18b is a massive hot Jupiter with M_P = 3.47 M_Jup and R_P = 1.31 R_Jup on a short-period orbit (P = 1.9 d) around a G9V star (Hébrard et al. 2011). The primary is likely to be a young star with an age of several hundred Ma, but the age is not well constrained. We find a secondary eclipse candidate near a phase of 0.47, corresponding to an eccentric orbit with θ_e = 0.1 ± 0.04. As with CoRoT-3b and CoRoT-13b, we obtain a high flux ratio corresponding to a median brightness temperature (2700 K). If we were to assume the eclipse to be real, the high eccentricity and temperature would both support the young stellar age hypothesis.

4. CoRoT-21b is a hot Jupiter with M_P = 2.26 M_Jup, R_P = 1.30 R_Jup on a P = 2.72 day orbit around a faint (V = 16) F8IV star (Pätzold et al. 2012). We find a marginal eclipse candidate with B₁₀ = 1.7. The signal passes the splitting test, but the mean brightness temperature estimate of 3100 K is again anomalously high. However, the flux ratio estimate has a large uncertainty, and the 95% credible intervals of T_br and T_eq overlap.

6.5. Shortcomings and possibilities for future improvement

While our brightness temperature estimates for CoRoT-1b and CoRoT-2b agree with the published values, all the new candidates show anomalously high brightness temperatures. For some of the candidates the temperatures could be explained by young age, atmospheric temperature inversion, deuterium burning, or tidal heating. However, the plausibility of different additional heat sources would need detailed modeling of the individual systems that is out of the scope of the current generalizing study aimed at presenting the new candidates.

Assuming normally distributed i.i.d. errors is a simplification, especially since the light curves are known to contain sources of systematic noise. The assumption of normality was tested to be a good approximation for the analyzed light curves, but the errors cannot be considered independent. Inclusion of correlated noise may increase our ability to distinguish true eclipse signals from the noise due to the stellar granulation, and will be considered for the future studies.

7. Conclusions

We have presented and characterized three new statistically significant eclipse candidates and four new marginal candidates, and offered an independent confirmation of the CoRoT-1b and CoRoT-2b eclipses.

All three of the new significant candidates have higher Bayes factors than the two confirmed eclipses, but they either imply properties that deviate from the ones derived from transits (e.g., the impact parameter of CoRoT-11b), or yield high brightness temperatures. Three of the four marginal candidates have higher Bayes factors than the CoRoT-2b eclipse. These candidates correspond to planets orbiting faint host stars and the light curves contain significant amounts of correlated noise, faintness and noise both contributing to the large uncertainty in the derived parameters.

The estimated optical bandpass brightness temperatures are higher than the expected equilibrium temperatures for all of our new candidates. This agrees with the results by Coughlin & López-Morales (2012), who find significant amounts of excess flux for their Kepler eclipses, and can be explained simply by the fact that planets are not black bodies. However, further observations of eclipse events in the infrared are necessary to confirm both the candidates and the characterization presented in this work.

Acknowledgments

The authors kindly thank the anonymous referee for the constructive and insightful comments. The present study was made possible thanks to observations obtained with CoRoT, a space project operated by the French Space Agency, CNES, with participation of the Science Program of ESA, ESTEC/RSSD, Austria, Belgium, Brazil, Germany and Spain. H.P. has received support from RoPACS during this research, a Marie Curie Initial Training Network funded by the European Commissions Seventh Framework Programme. H.P. and H.D. acknowledge funding by grant AYA2010-20982 of the Spanish Ministry of Economy and Competitiveness (MINECO). The author thankfully acknowledges the technical expertise and assistance provided by the Spanish Supercomputing Network (Red Espanola de Supercomputacion), as well as the computer resources used: the LaPalma Supercomputer, located at the Instituto de Astrofisica de Canarias.

References

All Tables

All Figures

	Fig. 1 Main results for CoRoT-1b. Shown are the eclipse model probability and the Bayes factor (a) and b), respectively); the lightcurve and the model, both folded over the period and binned with three different bin widths c); and posterior distributions for the flux ratio, eclipse center phase, eccentricity and the argument of the periastron (d)–g)).
In the text

	Fig. 2 As in Fig. 1, for CoRoT-2b.
In the text

	Fig. 3 As in Fig. 1, for CoRoT-6b.
In the text

	Fig. 4 As in Fig. 1, for CoRoT-11b.
In the text

	Fig. 5 As in Fig. 1, for CoRoT-15b.
In the text