Issue 
A&A
Volume 587, March 2016



Article Number  A149  
Number of page(s)  29  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201526297  
Published online  04 March 2016 
Inferring heat recirculation and albedo for exoplanetary atmospheres: Comparing optical phase curves and secondary eclipse data
^{1} Univ. Bordeaux, LAB, UMR 5804, 33270 Floirac, France
email: vonparisp@gmail.com
^{2} CNRS, LAB, UMR 5804, 33270 Floirac, France
Received: 10 April 2015
Accepted: 14 December 2015
Context. Basic atmospheric properties, such as albedo and heat redistribution between day and nightsides, have been inferred for a number of planets using observations of secondary eclipses and thermal phase curves. Optical phase curves have not yet been used to constrain these atmospheric properties consistently.
Aims. We model previously published phase curves of CoRoT1b, TrES2b, and HATP7b, and infer albedos and recirculation efficiencies. These are then compared to previous estimates based on secondary eclipse data.
Methods. We use a physically consistent model to construct optical phase curves. This model takes Lambertian reflection, thermal emission, ellipsoidal variations, and Doppler boosting, into account.
Results. CoRoT1b shows a nonnegligible scattering albedo (0.11 < A_{S} < 0.3 at 95% confidence) as well as small daynight temperature contrasts, which are indicative of moderate to high redistribution of energy between dayside and nightside. These values are contrary to previous secondary eclipse and phase curve analyses. In the case of HATP7b, model results suggest a relatively high scattering albedo (A_{S} ≈ 0.3). This confirms previous phase curve analysis; however, it is in slight contradiction to values inferred from secondary eclipse data. For TrES2b, both approaches yield very similar estimates of albedo and heat recirculation. Discrepancies between recirculation and albedo values as inferred from secondary eclipse and optical phase curve analyses might be interpreted as a hint that optical and IR observations probe different atmospheric layers, hence temperatures.
Key words: planets and satellites: individual: CoRoT1b / techniques: photometric / planets and satellites: atmospheres / planets and satellites: individual: TrES2b / planets and satellites: individual: HATP7b
© ESO, 2016
1. Introduction
In recent years, exoplanetary science has developed from a detectioncentered astronomical science towards a characterizationcentered planetary science. For basic properties of exoplanetary atmospheres such as albedo and heat redistribution from dayside to nightside, observational constraints are now available for a growing number of planets.
Fig. 1 Basic sketch of the geometry of the model. 
Radiative transfer and atmospheric modeling by, e.g., Sudarsky et al. (2000) predicts very low optical albedos for cloudfree hot Jupiters because of strong absorption by alkali metals. When silicate clouds form high enough in the atmosphere, however, optical albedos could be significantly higher (Sudarsky et al. 2000). The low optical albedos measured for a few hot Jupiters seem to confirm this (e.g., Rowe et al. 2006, 2008), whereas high albedos inferred for other planets seem to indicate the potential presence of clouds (e.g., Quintana et al. 2013; Demory et al. 2013; Esteves et al. 2013).
Theoretically, the most basic circulation patterns of hot Jupiters are relatively well understood (e.g., Showman & Guillot 2002; Madhusudhan et al. 2014; Heng & Showman 2015). Tidally locked planets in closein orbits always present the same hemisphere to the host star. In atmospheric models, the resulting strong irradiation contrasts and slow planetary rotations lead to inefficient recirculation throughout the IR photosphere of the planet (around pressures of 0.1–10 mbar). The incident stellar energy is reradiated before being circulated to the nightside. Consequently, strong temperature contrasts between dayside and nightside develop (e.g., Showman & Guillot 2002; Parmentier et al. 2013). When orbital distance or rotation period increases, atmospheric modeling predicts a transition to a circulation regime with much stronger recirculation and less pronounced daynight temperature differences (e.g., Showman et al. 2015). Furthermore, in most 3D atmosphere models of hot Jupiters, a strong equatorial zonal jet appears, that then results in a displacement of the hottest point away from the substellar point (e.g., Showman & Polvani 2011; PerezBecker & Showman 2013; Showman et al. 2015).
Observed thermal phase curves in the (near)infrared (IR) of a few hot Jupiters seem to confirm these predictions (e.g., Knutson et al. 2007, 2009; Stevenson et al. 2014). Assembling many different observations for a large number of hot Jupiters, Cowan & Agol (2011) and Schwartz & Cowan (2015) point out a possibly emerging trend, in line with theoretical predictions (e.g., PerezBecker & Showman 2013). Strongly irradiated planets show very inefficient recirculation, whereas for less irradiated planets, the whole range of recirculation is possible.
Cowan & Agol (2011) and Schwartz & Cowan (2015) use published secondary eclipse data (i.e., an estimate of the planetary dayside flux) and thermal phase curves to perform their homogeneous analysis. For a few planets, optical phase curves are available, obtained with the CoRoT and Kepler satellites. These are not taken into account in Schwartz & Cowan (2015). However, for hot planets, even optical phase curves offer some constraints on thermal radiation, and therefore albedo and heat recirculation.
So far, published studies of optical phase curves could not be used for such an analysis. This is, in each case, because of one of the three following reasons.
Firstly, some models do not take thermal emission into account (e.g., Mazeh & Faigler 2010; Barclay et al. 2012; Esteves et al. 2013, 2015). In such cases, attributing an observed phase curve asymmetry to an offset of the thermal hotspot is physically inconsistent (e.g., Esteves et al. 2015). Secondly, a few models (e.g., Snellen et al. 2009) only treat thermal emission and do not include scattered light. Therefore, inferring constraints on the scattering properties of the planet (as done in Snellen et al. 2009) is equally inconsistent. Thirdly, models that take both thermal and scattering components into account (e.g., Mislis et al. 2012; Faigler et al. 2013; Placek et al. 2014; Faigler & Mazeh 2015), do not use appropriate phase functions for these components. They assume, for instance, that thermal and scattered light have identical phase functions which is incoherent when assuming explicitly Lambertian scattering and blackbody thermal radiation (see below, Appendix A and Eqs. (7) and (17)).
Therefore, we present a new model with a physically consistent treatment of both components. We apply our model to three hot Jupiters with wellcharacterized IR and optical measurements, namely CoRoT1b, TrES2b and HATP7b, to compare to results from secondary eclipse analysis.
In Sects. 2 and 3, we describe the physical forward model, the inverse model and the Markov chain Monte Carlo (MCMC) approach used in this work. Section 4 describes the model setup and planetary scenarios, Sect. 5 presents the results, Sect. 6 a discussion, and we conclude with Sect. 7. The Appendices contain the fitting results, MCMC output, a short model verification as well as a discussion about Rayleigh scattering, nonLambertian phase functions and stellar parameter uncertainties.
2. Forward model
2.1. Geometry
The basic geometry setup of the starplanet system is shown in Fig. 1. We adopt a coordinate system where the substellar point is fixed throughout the calculations at 0° latitude and 180° longitude. Therefore, by definition, local latitude ϑ is 0° at the equator (range −90°<ϑ< 90°) and local longitude ϕ is 0° at local midnight (range 0°<ϕ< 360°). For planets with zero obliquity and a planetary rotation synchronized with the orbital period, our choice of coordinate system would establish a stationary map of the planet. For closein giant planets, this is a reasonable assumption (but see, e.g., Arras & Socrates 2010 or Rauscher & Kempton 2014 for a discussion).
The planetary “surface” is divided into cells of 2.5° × 2.5° size (72 cells in latitude, 144 in longitude). We define here “surface” as the optical photosphere, i.e., where the planetary atmosphere becomes optically thick to visible radiation. It is generally identified with the surface of the sphere with radius R_{P}. The number of cells is reasonable for computational purposes and still allows for smooth light curves without noticeable effects of discretization. For each cell, the local surface normal n, stellar and observer directions (s and o, respectively) are calculated. A cell contributes to the reflected light if both (i.e., dayside) and (i.e., visible to the observer) are satisfied. The nightside is defined with , and correspondingly, the part of the nightside visible by the observer with and , simultaneously.
The subobserver latitude θ is given by the inclination i of the orbital plane with respect to the observer: (1)i.e., an edgeon orbit has i = 90°, and a faceon orbit has i = 0°. The subobserver longitude φ as a function of time t is given by the orbital phase α (α = 0 and t = 0 at primary transit, or equivalently, inferior conjunction, see Fig. 1): (2)with (3)where T is the true anomaly and ω_{P} the argument of periastron. The true anomaly is calculated from the eccentric anomaly E: (4)with e eccentricity. E is determined by numerically solving Kepler’s equation: (5)with M mean anomaly. The solution is obtained with a publicly available Fortran procedure^{1}, based on a method described in Meeus (1991). The mean anomaly M is given by (6)where , P_{orb} is the planetary orbital period, t_{peri} is the time of periastron passage and ⌊ x ⌋ represents the floor function, i.e., the greatest integer less than or equal to x^{2}.
2.2. Reflected light
Stellar light incident on the planetary dayside is partly reflected back into space and towards the observer. The amount of reflected light reaching the observer depends, for instance, on the scattering properties of the planet (e.g., Rayleigh scattering produces a different phase function than Mie scattering, see, e.g., Madhusudhan & Burrows 2012 for a review). In this work, in the absence of any reliable information on the scattering properties, the Lambert approximation of diffuse scattering is used. In Appendix C we discuss the possible influence of Rayleigh scattering and other phase functions on the phase curve.
In the Lambert approximation, for each cell contributing to the observed flux (), the flux (in W m^{2}) received by the observer from this cell, F_{r,o,cell}, is given by (7)with z_{s} the local stellar zenith angle, z_{o} the local observer zenith angle, F_{∗ ,p} the stellar flux at the planet’s orbit (in W m^{2}), ΔS the surface element of the cell on the planet (in sr m^{2}), d the observerplanet distance and A_{S} the (potentially wavelengthdependent) planetary scattering albedo. A_{S} is assumed to be constant with time, hence neglecting any timedependent processes such as cloud formation etc. Note that, since we use the Lambertian approximation in Eq. (7), the scattering albedo is related to the geometric albedo A_{G} in a simple manner: (8)The surface element ΔS of the cell, as seen from the planet’s center, is calculated as (9)with R_{P} the planetary radius and ΔΩ the solid angle (in sr) of the cell (angular extent 2.5° × 2.5°, see above).
The stellar flux at the planet’s orbit, F_{∗ ,p}, is calculated as follows: (10)where I_{∗ ,s} is a stellar model intensity (in W m^{2} sr^{1} μm^{1}) at the stellar surface, q_{I} is the instrumental filter function, λ_{low} and λ_{high} define the wavelength interval of the bandpass, r the starplanet distance and R_{∗} is the stellar radius^{3}. The numerical integration for Eq. (10) is done with a standard trapezoidal integration scheme using a roughly 1 nm spectral resolution.
Stellar intensities I_{∗ ,s} are obtained from the ATLAS stellar atmosphere grid (Castelli & Kurucz 2004)^{4}. Instrumental filter functions T_{l} are taken, for CoRoT1, from Snellen et al. (2009) and, for TrES2 and HATP7, the Kepler Handbook^{5}.
The timedependent planetstar distance r is calculated with (11)with a the semimajor axis.
The total reflected flux F_{r,o} received by the observer is the sum over all contributing cells: (12)
2.3. Emitted light
In addition to reflected starlight, the planet also emits thermal radiation that can contribute to the overall phase curve. In optical phase curves, this contribution would be negligible for longperiod (and consequently colder) planets. However, for closein hot Jupiters, the thermal component can become comparable to (or even dominate) the reflected component of the phase curve.
In this work, we make the specific assumption of two hemispheres (nightside and dayside) which radiate as a blackbody with uniform temperatures T_{night} and T_{day}, respectively.
A widely used approach (for example, Snellen et al. 2009; Alonso et al. 2009; Mislis et al. 2012; Schwartz & Cowan 2015) is to relate these temperatures to the Bond albedo A_{B} and the efficiency of heat redistribution ϵ. Temperatures are then calculated, based on Cowan & Agol (2011): (13)and (14)where T_{∗} is the stellar effective temperature and the additional factor (1−e^{2})^{0.5} accounts for the mean flux received over an orbit (Williams & Pollard 2002).
The parameter ϵ is a reparameterization of the geometrical redistribution factor f. Spiegel & Burrows (2010) define f via the apparent dayside flux F_{d} (i.e., close to secondary eclipse) and the total planetary luminosity L_{P}: (15)For perfect heat recirculation, the entire planet is at a uniform temperature, and thus , because both dayside and nightside hemispheres contribute to the planetary luminosity. At zero heat recirculation, the nightside emission is zero, and each point on the dayside hemisphere emits with its local radiative equilibrium temperature. Hence, as shown by Spiegel & Burrows (2010) and Cowan & Agol (2011), .
Assuming radiative equilibrium, meaning that the total stellar flux F_{S} intercepted by the planet equals its luminosity L_{P}, and approximating the star as a blackbody, one can rearrange Eq. (15) to yield (16)For ϵ to vary between 0 and 1 (corresponding to no recirculation and perfect redistribution, respectively) then simply requires a linear transformation of f, which leads directly to Eq. (13).
Physically, ϵ is determined by the atmospheric circulation and the strength of winds and jets that transport heat away from the illuminated hemisphere. For strongly irradiated planets, the radiative timescale is expected to be much shorter than the advective (dynamical) timescale, therefore large daynight temperature contrasts and small values of ϵ are expected. For less irradiated planets, a range of circulation regimes is possible (e.g., Showman & Guillot 2002; Showman & Polvani 2011; PerezBecker & Showman 2013; Heng & Showman 2015; Showman et al. 2015).
Another option in the model is to retain T_{night} and T_{day} as free parameters for the inverse modeling, an approach used by, e.g., Placek et al. (2014).
Since blackbody radiation is isotropic, the thermal flux F_{t,o,cell} received by the observer from each cell is given by (17)with B(T_{c},λ) the blackbody intensity at the cell’s temperature T_{c} and T_{c} = T_{night} or T_{c} = T_{day}, depending on the location of the cell.
Again, the total emitted flux F_{t,o} is the sum over all cells which are visible for the observer (i.e., ): (18)Note that in our approach, thermal emission produces a different phase curve behavior compared to Lambert scattering of stellar radiation because of the additional factor cosz in Eq. (7) compared to Eq. (17). Therefore it is potentially possible to disentangle reflected from emitted light. This effect is illustrated in Fig. 2.
Fig. 2 Effect of thermal emission on the phase curve. See text for discussion. 
Figure 2 shows the reflective and the thermal contribution for a hot Jupiter planet in a 2day orbit around a Sunlike star (500–1000 nm bandpass, A_{B} = A_{S} = 0.15, ϵ = 0). In this example, the planetary mass is set to zero for illustration purposes. The reflected component is shown as a dashed line, the thermal component of the phase curve is the dotted line. Also shown are the combined phase curve (plain line) and a purely reflective phase curve (dotdashed line) that has the same maximum amplitude as the combined phase curve. Due to the additional angle dependence of the reflected light (cosz_{s} in Eq. (7)), the reflected phase curve is steeper than the thermal one and shows a more pronounced peak towards phase α = π. Based on the slope, if the orbital period is short enough and the photometric precision good enough, it is possible to determine dayside and nightside temperatures, hence Bond albedo and heat redistribution, contrary to what is generally assumed in previous studies (e.g., Esteves et al. 2013, 2015; Schwartz & Cowan 2015) that did not incorporate thermal radiation consistently.
Of course, note that in reality, nonuniform distribution of temperatures (due to, e.g., chemical composition changes, Agúndez et al. 2012) will complicate the interpretation of thermal phase curves. However, because of the currently limited data, assuming uniform hemispheres is justifiable.
2.4. Ellipsoidal variations and Doppler boosting
In addition to the planetary contributions to the phase curve, our model also considers two modulations of the stellar light induced by the planet. These are the ellipsoidal variations (e.g., Pfahl et al. 2008) and the Doppler boosting (e.g., Loeb & Gaudi 2003). In short, ellipsoidal variations are due to the tidal deformation of the star by the orbiting planet. As a consequence, the star represents a varying cross section to the observer, hence, the luminosity changes periodically, with a period half that of the orbital period of the companion. Doppler boosting is a consequence of the Doppler shift of the stellar spectrum induced by the stellar reflex motion. As the star orbits the barycentre of the starplanet system, the emitted stellar light shifts its wavelength, and thus periodically more or less light is emitted in the bandpass of the used instruments.
We use the formalism developed in Quintana et al. (2013) to calculate the respective contrasts:
where A_{ell}, A_{dopp} are the amplitudes of the ellipsoidal and Doppler variations. Note that we do not fit for a potential phase lag between tidal bulge on the star and the planet (as in, e.g., Barclay et al. 2012). Furthermore, we do not incorporate other harmonics of the orbital period into our ellipsoidal variation term (contrary to what was as done, e.g., in Esteves et al. 2013). A_{ell} is given by (see also Eq. (2) in Quintana et al. 2013): (21)with M_{∗}, M_{p} the stellar and planetary mass, respectively. α_{ell} is a parameter determined by the stellar limb (u) and gravity (g) darkening: (22)The coefficients u and g depend on stellar characteristics such as effective temperature T_{∗}, metallicity [ Fe / H ] and surface gravity log g_{s} (determined from M_{∗} and R_{∗}). As in Barclay et al. (2012), Quintana et al. (2013) or Esteves et al. (2013), we use precalculated tables for u and g to interpolate linearly in T_{∗}, [ Fe / H ] and log g_{s}. These tables are taken from model calculations presented in Claret & Bloemen (2011) for the Kepler and CoRoT bandpasses. We adopt their coefficients obtained with a microturbulence velocity of 2 km s^{1} with ATLAS stellar models and, in the case of u, fitted with a linear least squares approach.
A_{dopp}is given by (see also Eq. (5) in Quintana et al. 2013): (23)where c is the speed of light, K the radial velocity semiamplitude and α_{dopp} a parameter which depends on the wavelength of observation λ_{obs} and the stellar effective temperature. As in Quintana et al. (2013), we use the approximate equations from Loeb & Gaudi (2003) for α_{dopp}: (24)where (h Planck’s constant, k_{B} Boltzmann’s constant). This approach of calculating (3−α_{dopp}) is somewhat different from the approach used in, e.g., Esteves et al. (2013, 2015). However, comparing results for TrES2b (their work, 3.71, to 3.87 with our approach) and HATP7b (3.41 to 3.59) indicates that both approaches yield values that are within 5%. Hence mass determinations are expected to be comparable.
Kis determined as (with G the gravitational constant) (25)
2.5. Asymmetries
As demonstrated by, e.g., Demory et al. (2013) and Esteves et al. (2015), many exoplanets show asymmetries in their phase curves, with respect to secondary eclipse. This means that the maximum amplitude is reached before (after) secondary eclipse, meaning that the maximum brightness is shifted eastwards (westwards) of the substellar spot. In the case of a westwards shift, it has been interpreted as the presence of clouds that enhance the reflectivity of the “morning”^{6} side. When the shift is eastwards, i.e., the “evening” side is brighter, previous studies attributed this to a shift in the hottest region of the atmosphere. Such a shift is associated with atmospheric circulation which transports heat away from the substellar point before it is reradiated. Most 3D atmospheric models of hot Jupiters produce such an offset in thermal emission (e.g., Showman & Guillot 2002; Showman et al. 2015), and IR phase curves seem to confirm the theoretical predictions (e.g., Knutson et al. 2007).
This change in brightness can be accounted for in the model in two ways. First, for the reflectedlight component of the phase curve, we implement a simple darkbright model, similar to the approach chosen in Demory et al. (2013). Part of the planet has a scattering albedo A_{S} (see Eq. (7)), and part of the planet has a different scattering albedo d_{S}·A_{S}, with d_{S} being a free parameter such that d_{S}·A_{S} ≤ 1. The extent in longitude of the latter part of the planet is controlled by two further parameters, namely l_{start} and l_{end}. We do not consider any latitudinal variation of the scattering properties.
Fig. 3 Illustration of the asymmetric models. Left: reflective dayside (red) with bright “morning” (green). Right: thermal offset modeled as a shifted dayside (red). See text for discussion. 
Second, for the thermal component of the phase curve, we consider a simple offset Θ_{d} of the dayside such that the dayside has an extent in longitude between Θ_{d} and 180° + Θ_{d}. Figure 3 illustrates the asymmetric models.
2.6. Phase curve
The final timedependent contrast C between star and planet is then (26)The stellar flux at the observer, F_{∗ ,o}, is calculated analogous to F_{∗ ,p} (see Eq. (10)): (27)
2.7. Model parameter summary
In total, our full physical model as described above contains up to 19 parameters, listed in Table 1.
Parameters of the forward model.
These 19 parameters however are not independent. For instance, if A_{B} and ϵ are fixed, T_{d} and T_{n} can be determined using Eqs. (13) and (14).
3. Inverse model
Planetary scenarios for the forward model.
We use the Bayesian formalism to calculate posterior probability values p(V_{P}  D) for the parameter vector V_{P} in the model, given a set D of observations. (28)The likelihood p(D  V_{P}) is calculated assuming independent measurements and Gaussian errors for the individual data points. The priors p(V_{P}) are taken to be uninformative over the entire parameter range allowed (for example, uniform over [0, 1] for albedo and heat redistribution).
3.1. MCMC algorithm
To sample the full parameter space, we adopt a MCMC approach. In this work, we use the emcee python package developed by ForemanMackey et al. (2013), implementing an algorithm described in Goodman & Weare (2010). emcee uses multiple chains (in this work, 500–1000) to sample the parameter space. The algorithm proposes, for each chain, new positions based on the position of the entire ensemble of chains. Compared to more traditional MCMC approaches, emcee converges quicker and is less likely to be dependent on initial conditions. Also, when using a high number of chains, the algorithm is less likely to become stuck in local minima since it is possible to eliminate chains from the ensemble (e.g., Hou et al. 2012).
To ensure good convergence and avoid any contamination by initial conditions, the chains were run for 500–2000 steps (>10–20 autocorrelation lengths for each parameter). The first few autocorrelation lengths were considered as burnin and discarded for the calculation of parameter uncertainties. Convergence was checked by inspecting visually the evolution of the mean of the entire ensemble and calculating the GelmanRubin test. Initial positions were obtained with a random sample within the assumed prior to allow the sampler to start by exploring the entire parameter space.
Fig. 4 Trace plots of model parameters for the HATP7b “standard + asy” model (see Table 2). Blue line traces the ensemble median, red lines correspond to the [0.16, 0.84] percentiles (the dark red region in Fig. 5 below). 
Uncertainty ranges are calculated by marginalizing over the posterior distribution, thinned by 60 steps, covering roughly 1–3 autocorrelation lengths of the particular parameter in question. We then determine 68% and 95% credibility regions as the [0.16, 0.84] and [0.03, 0.97] mediancentered percentiles, respectively, of the cumulative probability distributions (CDF). If the parameter distribution were to be Gaussian, these credibility regions would correspond to the 1 and 2σ uncertainties, respectively. Figure 5 illustrates our method to determine credibility regions and bestfit parameters. Bestfit parameters are determined as the maximum a posteriori (MAP) set of parameters, i.e., the walker with the highest a posteriori probability at the last step of the algorithm. Note that the MAP does not correspond to the median of the CDF in most cases. Furthermore, depending on the nature of the likelihood surface sampled by the chains (degeneracies, slope, etc.), the MAP, as determined from our sample, can occasionally even lie outside the [0.03, 0.97] percentile (see Tables D.1–D.3, Figs. D.1–D.8).
Fig. 5 Example of determination of uncertainty ranges via the cumulative probability distribution: Scattering albedo of CoRoT1b in the “standard” scenario (see text for further details). 68% credibility region in dark red, 95% credibility region in light red. Dashed lines indicate maximum a posteriori (MAP) value. 
3.2. Goodnessoffit criteria and model comparison
We use three standard quantities (e.g., Feigelson & Babu 2012) to evaluate the bestfitting models obtained from the MCMC model: χ^{2}, the reduced and the Bayesian Information Criterion (BIC). These are evaluated for the MAP parameter set (not the median of the CDF).
χ^{2}, the sum of the weighted squared residuals, is defined as (29)where N_{D} is the number of data points, σ the corresponding uncertainties and M_{VP} is the model prediction for the parameter vector V_{P}.
The reduced is generally calculated by (30)where N_{P} is the number of parameters. In order to be considered a good fit, should be of the order of unity.
Note that by adding more parameters to the fitting model, the χ^{2} will decrease. Therefore, we use another criterion, the BIC (valid when N_{D} ≫ N_{P}, which is the case for our calculations), defined as (31)The model that minimizes the BIC is taken to be the preferred model. The BIC penalizes overly complex models when complexity or sophistication is not warranted by the data, and also allows direct model comparison. For instance, the model probability ratio p_{M} (i.e., the probability that the model is preferred over another one), can be expressed as (32)where ΔBIC = BIC_{M1}−BIC_{M0} is the difference between the model M_{1} under consideration, and the model M_{0} that minimizes the BIC.
The BIC is an approximation to the evidence (or marginal likelihood) E_{model} of a given model: (33)The calculation of the integral in Eq. (33) is in most case computationally very expensive. In emcee, this integral is estimated using a socalled thermodynamic integration, based on an algorithm proposed by Goggans & Chi (2004). However, the calculation is very timeconsuming (up to ~50 times longer than the MCMC sampling), and does not, in our cases, provide any qualitatively additional information, compared to the BIC. Therefore, we only report the BIC and base our discussions on Eq. (32).
4. Model setup
4.1. Planets
4.1.1. CoRoT1b
CoRoT1b is the first planet discovered from space (Barge et al. 2008) and the first planet with a detected optical phase curve (Snellen et al. 2009). Secondary eclipses of CoRoT1b have been detected in the optical (Alonso et al. 2009) and IR, both from ground (e.g., Rogers et al. 2009; Gillon et al. 2009) and from space (e.g., Deming et al. 2011). These studies suggested that CoRoT1b would be a very dark planet, with a low albedo (A_{G}≲0.1) and inefficient heat redistribution (ϵ ≲ 0.2), leading to high dayside temperatures of the order of 2300 K.
We use the binned CoRoT redchannel phase curve data presented by Snellen et al. (2009, Fig. 1). All model parameters (stellar parameters, orbital period and inclination, planetary mass and radius) are taken from Barge et al. (2008), similarly to what was done in Snellen et al. (2009). The orbit is assumed to be circular. This is consistent with the analysis of secondaryeclipse timing (e.g., Rogers et al. 2009; Alonso et al. 2009; Deming et al. 2011) that constrain the eccentricity to values of e ≲ 0.03. Given the data quality of the optical phase curve, the results are not sensitive to eccentricity. When allowing e and ω_{P} to be fitted (simulations that are not shown here), constraints on either A_{B} or ϵ did not change appreciably compared to the circular case.
4.1.2. TrES2b
TrES2b is a transiting hot Jupiter (O’Donovan et al. 2006) discovered by the TrES survey. It was the first transiting exoplanet discovered in the Kepler field prior to Kepler’s 2009 launch. Groundbased and Spitzer secondary eclipse photometry suggests moderately high temperatures around 1500 K (e.g., O’Donovan et al. 2010; Croll et al. 2010). The optical phase curve has been detected in the Kepler data (e.g., Kipping & Spiegel 2011; Barclay et al. 2012; Esteves et al. 2013). Results indicate a very low geometric albedo (A_{G} ≤ 0.02) and a rather efficient daynight redistribution of energy (ϵ>0.5), since dayside and nightside temperatures are quite similar (around 1300–1500 K, Esteves et al. 2013).
We use the binned Kepler phase curve data (Barclay et al. 2012, their Fig. 5). All model parameters (stellar parameters, orbital period and inclination, planetary radius) are taken from Barclay et al. (2012). The orbit is assumed to be circular, since secondaryeclipse timing provided constraints consistent with zero eccentricity (e.g., O’Donovan et al. 2010; Croll et al. 2010) and optical phase curve analysis by Esteves et al. (2015) did not find any significant eccentricity.
4.1.3. HATP7b
As TrES2b, HATP7b is a very hot Jupiter discovered in the Kepler field (Pál et al. 2008) prior to the launch of the satellite. It was the first Kepler planet with a measured optical phase curve (Borucki et al. 2009). Further analysis of the phase curve demonstrated a detection of ellipsoidal variations (Welsh et al. 2010). Christiansen et al. (2010) used Spitzer secondary eclipse data to infer maximum brightness temperatures of more than 3000 K. Esteves et al. (2013) presented a new optical phase curve using full 3year Kepler photometry. Measured phase curve amplitudes and secondary eclipse depths vary considerably between Borucki et al. (2009), Welsh et al. (2010) and Esteves et al. (2013), leading to differences in inferred dayside and nightside temperatures of the order of 500 K which are most likely a consequence of the extended data set in Esteves et al. (2013). Similar brightness temperature differences of 300–500 K for the Spitzer secondary eclipse data have been found by Cowan & Agol (2011), compared to Christiansen et al. (2010), which they attribute to the use of different stellar models. Based on calculated brightness temperatures and optical phase curves, Christiansen et al. (2010) inferred a very inefficient heat redistribution (ϵ close to zero) between dayside and nightside and modest geometric albedos (A_{G}< 0.1). Esteves et al. (2013) found a geometric albedo of A_{G} = 0.18 and a relatively homogenous temperature distribution, in contrast to Christiansen et al. (2010). Schwartz & Cowan (2015) found moderate albedos, slightly lower than Esteves et al. (2013) and moderate recirculation values, also in contrast to previous analysis (Christiansen et al. 2010).
We use the binned Kepler phase curve data (Esteves et al. 2013, their Fig. 3). However, model parameters are not taken from Esteves et al. (2013) or Esteves et al. (2015). Instead, stellar parameters are taken from Van Eylen et al. (2012), and planetary parameters (inclination, radius) are taken from Van Eylen et al. (2013). We refer to Appendix B for a discussion on our choice of parameters. Again, we assume a circular orbit, since detailed radialvelocity data is consistent with e = 0 (e.g., Pál et al. 2008; Winn et al. 2009). Also, previous phase curve analysis did not find a hint of significant eccentricity (Esteves et al. 2015).
4.2. Photometric fits
Snellen et al. (2009) analyze the phase curve of CoRoT1b in terms of the normalized flux F_{norm}, normalized to the primary, instead of secondary, eclipse (i.e., without transit, the flux at zerophase would be unity). Therefore, we calculate the forward model as: (34)Both Barclay et al. (2012) and Esteves et al. (2013) fit the phase curves in terms of variations in the photometric light curve, as in Eq. (26). They however introduce another, nonphysical parameter, a zeropoint flux offset f_{0}. This parameter is related to the data reduction and not a priori linked to any physical characteristics of the starplanet system. The light curve L_{C} is then described with the following equation: (35)We will use this equation for our analysis of TrES2b and HATP7b. Hence, f_{0} is added to the tally of free parameters of the physical model (Table 2).
4.3. Planetary scenarios
As summarized in Table 2, we explore several different scenarios for the three planets. The main difference between these models lies in the treatment of the thermal component of the phase curve.
The first scenario (“standard”) assumes scattering albedo A_{S} and heat redistribution as fitting parameters and uses Eqs. (13) and (14) to calculate hemispheric temperatures. The Bond albedo A_{B} is fixed to the scattering albedo (A_{S} = A_{B}). This allows our results to be compared directly to previous inferences of albedo and heat recirculation (e.g., Snellen et al. 2009; Schwartz & Cowan 2015) and is consistent with previous studies (e.g., Alonso et al. 2009). The equality between Bond albedo and scattering albedo is motivated by the fact that a significant fraction a_{V} of stellar light is emitted in the CoRoT red channel or the Kepler bandpass (a_{V} ≈ 0.3 for CoRoT1).
The second scenario (“free albedo”) relaxes this tight coupling between Bond albedo and scattering albedo. We allow both A_{B} and A_{S} to vary freely, but still calculate dayside and nightside temperatures from the Bond albedo, as in the first scenario. Based on energy conservation, however, we put some constraints on the Bond albedo: (36)This equation takes into account the contribution of the scattering albedo to the overall radiative budget. The lower limit of the Bond albedo (A_{B,low} = a_{V}·A_{S}) is a hard lower limit since it implies zero albedo outside the bandpass. Similarly, when assuming an albedo of unity outside the bandpass, we obtain the strict upper limit of A_{B,high} = a_{V}·A_{S} + (1−a_{V}).
In a third approach (“free T”), we use T_{d} and T_{n} as fitting parameters (instead of Bond albedo and heat recirculation) and only impose T_{n} ≤ T_{d}.
To compare directly to the phasecurve analysis by Snellen et al. (2009), we then use a fourth model approach for CoRoT1b. We set A_{S} = 0, i.e., the phase curve is produced by thermal emission only (“no scattering”). This was done because the physical model used in Snellen et al. (2009) only accounts for thermal radiation (their Eq. (1)).
All scenarios for CoRoT1b assume symmetric phase curves and do not fit for planetary mass because of the relatively low signaltonoise ratio and reduced phase resolution of the binned phase curve.
For TrES2b and HATP7b, however, the data quality is good enough that ellipsoidal variations are clearly seen. Therefore, we take the planetary mass to be a fit parameter. We also fitted the data with asymmetric standard models (see Table 2), either with asymmetric scattering, a dayside offset or a combination of both.
4.4. Energy balance
The energy balance can also be expressed in terms of received and emitted flux. For this, we divide the planet in uniform hemispheres (see Sect. 2). The Bond albedo determines the overall received flux which must then be reemitted on the day and night hemispheres. (37)where the sums contain the measured fluxes in the various wavelength bands (N_{day} bands of the dayside spectrum, N_{night} of the nightside spectrum). Hence, under the assumption that outside the measured bands, no flux is emitted, we obtain an upper limit for the Bond albedo. This constraint is physically sound and does not rely on specific assumptions except radiative equilibrium and the hemispheric uniformity.
5. Results
5.1. Convergence results
The results of the convergence tests and the trace plots for all parameters are shown in Appendix E. All simulations seem to have reached a stationary distribution and thus converged. Most parameters also pass the GelmanRubin (GR) test (generally, for most MCMC tools, a value of less than 1.1–1.2 is considered acceptable). Note, however, that the GR test is a test for good mixing of the ensemble. Thus, even when a stationary distribution is reached because of a large number of chains used in the calculation, parameters might fail the GR test (values larger than 1.2). This usually means that the interchain variance is large compared to the variance of individual chains (in our cases, sometimes orders of magnitude). For strongly nonlinearly correlated parameters (e.g., in the “standard + off” and “standard + both” scenarios of HATP7b, see below), it will take a long time for a single chain to explore the entire permitted range. Thus, the GR statistic will be large, without necessarily impacting the convergence of the ensemble. This is clearly shown by the fact that the MCMC algorithm actually finds strong correlations.
A particularly good example of this are the “standard + asy” and “standard + both” scenarios for TrES2b (see Figs. E.4 and E.5, Table E.2). Since the parameters describing the albedo asymmetry are tightly anticorrelated, the MCMC algorithm finds basically two separated solutions (see also Fig. D.5), resulting in “bad” GR values. When restricting the calculations to one of the solutions (not shown), the GR test is passed by all parameters (values less than 1.17 in the “both” scenario, less than 1.09 in the “asymmetric” scenario).
5.2. CoRoT1b
Figure 6 shows the data and the bestfit model from Snellen et al. (2009) as well as our bestfit models from the various planetary scenarios (Table 2). Since the “free T” and the “free albedo” bestfit models are virtually identical, only the latter is shown, for clarity. Clearly, the standard and freealbedo bestfit models do not differ by much. It is also apparent that both the models presented in this work and the model of Snellen et al. (2009) provide reasonable fits to the data.
Fig. 6 CoRoT1b redchannel phase curve: Comparison of bestfit models with data (red) and fit by Snellen et al. (2009; yellow). Primary transit and secondary eclipse not shown. 
Bestfit parameters as well as goodnessoffit criteria and MCMC posterior parameter distributions are reported in Appendix D (Table D.1, Figs. D.1 and D.2). Based on the BIC value, the data seem to slightly favor the standard scenario, however, all models in Table D.1 seem acceptable.
Furthermore, fit results suggest that the phase curve is dominated by scattering rather than by thermal emission. This can be seen in Table D.1 from the fact that the inferred value of the scattering albedo A_{S} is largely unaffected by the choice of the thermal model (0.11 <A_{S}< 0.3 at 95% confidence in the “free A” model). The combined arithmetic mean of the scattering albedo in the three models (“standard”, “free A”, “free T”) is A_{S} = 0.22. The independence of A_{S} of the thermal model is illustrated in Fig. 7.
Fig. 7 Marginalized posterior distributions for CoRoT1b scattering albedo in different models. A_{S} is approximately independent of the thermal model (0.11 <A_{S}< 0.3 at 95% confidence in the “free A” model). 
Constraints for Bond albedo are weak, and heat recirculation is essentially unconstrained. In addition, constraints on ϵ show some dependence on the choice of the thermal model employed (see Fig. 8). This suggests that indeed the optical phase curve does not contain much information on the thermal component, hence is dominated by reflected starlight rather than thermal emission.
Independent circumstantial evidence for this can be drawn from the observed, flat transmission spectrum of CoRoT1b (Schlawin et al. 2014) which can be interpreted as a hint for the presence of clouds (or, at least, a reflecting layer in the upper atmosphere).
Fig. 8 Marginalized posterior distributions for CoRoT1b ϵ (left) and A_{B} (right) in different models. Constraints on A_{B} are weak. ϵ is unconstrained and depends on the thermal model. 
Figure 9 shows the constraints on dayside and nightside temperatures in the “free T” model. Inferred dayside temperatures are much lower than the IR brightness temperatures derived from secondaryeclipse measurements. Essentially, results are consistent with zero phase curve contribution from thermal emission, in accordance with the “standard” and “free A” models. Also, in accordance with Snellen et al. (2009), we find that the nightside emission is consistent with zero.
Fig. 9 Marginalized posterior distributions for CoRoT1b dayside and nightside temperatures in the “free T” model. 
Fig. 10 Joint credibility regions of recirculation and geometric albedo (see Eq. (8)) for CoRoT1b in the “standard” (red dots) and “no scattering” (blue dots) scenarios. Orange contour: 1σ uncertainty region in Schwartz & Cowan (2015). Our “standard” model strongly disagrees with previous work. 
Our results are somewhat contrary to the results from the IR secondaryeclipse measurements (Schwartz & Cowan 2015) and conclusions of the analysis of the optical phase curve by Snellen et al. (2009). Both studies used the same formalism as our “standard” scenario (i.e., using Eqs. (13) and (14)), and they conclude that CoRoT1b is probably a lowalbedo planet with inefficient heat recirculation. Note, however, that the phase curve model used by Snellen et al. (2009) only takes thermal radiation into account, hence the scattering albedo is, by default, zero. Our results suggest significant scattering and, depending on the thermal model, at least some recirculation.
These contrasting results are illustrated in Fig. 10. For this, we calculated the joint credibility regions in the twodimensional A_{G}ϵ space such that points simultaneously lie in the 95% credibility regions of each parameter. The joint credibility region thus corresponds to an approximately 90% probability. Figure 10 shows these credibility regions for both the “standard” and the “no scattering” scenario. It is clear that the 1σ uncertainty region of Schwartz & Cowan (2015) and our joint credibility region barely overlap in the “standard” case. By contrast, the “no scattering” case yields approximately the same constraints as the analysis by Snellen et al. (2009) and Schwartz & Cowan (2015).
However, the noscattering case is equivalent to imposing a strong prior on the scattering albedo (A_{S} = 0) that seems rather adhoc. Therefore, on physical grounds, we prefer our standard model over the noscattering case, even though goodnessoffit criteria could not be used to decide formally which model to prefer, given that the ΔBIC is small (see Table D.1).
Despite the lack of information on the thermal emission of CoRoT1b from its optical phase curve, we can however still put some constraints on the Bond albedo using the estimated value of A_{S} and Eq. (36). Fit results for the freeT scenario (see Table D.1) translate, at 95% confidence, to 0.03 <A_{B}< 0.82, or, when using the bestfit values, 0.06 <A_{B}< 0.8, with a_{V} ≈ 0.26 which is consistent with results from the “free A” scenario.
When using the reported IR brightness temperatures of CoRoT1b (Rogers et al. 2009; Deming et al. 2011) and our optical brightness temperatures (N_{day} = 4, N_{night} = 1 in Eq. (37)), we obtain A_{B}< 0.85, consistent with results from Eq. (36). This is not a strong constraint, since the spectral coverage is not large. However, it is a mostly modelindependent result. Especially, nightside emission measurements are missing, except for our optical phase curve analysis (resulting in a nondetection since the nightside emission is consistent with zero at the level of the measurement errors).
5.3. TrES2b
Figure 11 shows the various bestfit models of the different MCMC scenarios. It is apparent that the main difference between symmetric and asymmetric models is in the second peak where asymmetric models provide a better fit to the data. For clarity, the “free A” and “free T” scenarios are not shown. In Table D.2, we state bestfit parameters as well as 95% credibility regions for the parameters. Parameter posterior distributions are shown in Figs. D.3–D.5.
Fig. 11 TrES2b phase curve: comparison of bestfit models with data (red) and fit by Barclay et al. (2012; orange). Primary transit and secondary eclipse not shown. 
Results presented in Fig. 12 confirm the very dark nature of TrES2b, with scattering albedos A_{S}< 0.03 at 95% credibility, as already inferred by previous authors (e.g., Kipping & Spiegel 2011; Barclay et al. 2012; Esteves et al. 2013, 2015).
Fig. 12 Marginalized posterior distributions for TrES2b A_{S} in different models. TrES2b is a dark planet in all models (A_{S}< 0.03 at 95% confidence). 
Furthermore, as shown in Fig. 13, the value of ϵ depends strongly on the choice of the planetary scenario (asymmetric vs. symmetric), as is the case for CoRoT1b. For the asymmetric models, this is immediately obvious. With ϵ close to unity, there will be no strong contrast between dayside and nightside, hence no asymmetry can be produced in the “standard + off” scenario. In order to produce a noticeable effect on the phase curve, inefficient heat recirculation is required (ϵ close to zero). As can be seen in Fig. 13, the distributions for ϵ in the “standard + off” and “standard + both” scenarios are close to each other, which suggests that the preferred mechanism to produce an asymmetric phase curve is probably thermal radiation, not reflected light (i.e., ϵ is probably small).
Fig. 13 Marginalized posterior distributions for TrES2b ϵ in different models. ϵ depends strongly on the chosen model. 
Figure 14 shows the marginalized posterior distributions for the inferred planetary mass. Also shown are results of previous phase curve modeling by Barclay et al. (2012) and Esteves et al. (2013) as well as RV mass determinations. It is obvious that the precision of the photometric mass is worse than that of the RV data. However, our mass values are consistent with previous studies.
Fig. 14 Marginalized posterior distributions for TrES2b M_{P} in different models. Mass from previous studies (including RV measurements) in gray. 
It is clear that our model provides a relatively good fit to the observed phase curve (Fig. 11, –2.2, see Table D.2). Compared to the bestfit model by Barclay et al. (2012), our symmetric models consistently calculate a noticeably higher photometric contrast posteclipse. However, as already noted by Esteves et al. (2015), this is simply because the model of Barclay et al. (2012) allows for a separate, independent fitting of the beaming and ellipsoidal amplitudes that are adjusted to compensate for the apparent decrease in the phase curve (see also discussion in Faigler & Mazeh 2015). This is also the main reason why beaming and ellipsoidal masses do not agree with each other in Barclay et al. (2012) or Esteves et al. (2013).
When allowing for asymmetric phase curves, the fit becomes slightly better. In terms of the respective BIC values (see Table D.2), the “standard + off” scenario is slightly favored, although a ΔBIC ≈ 3.5 is not enough to detect firmly an asymmetry. Our tentative detection is therefore not in contradiction with conclusions of Esteves et al. (2015) who state that symmetric models are favored for TrES2b.
Figure D.4 (right panel) shows some interesting correlations between ϵ, A_{S} and Θ_{D}. For an increasing albedo, the offset of the dayside also increases. This is because, with increasing contribution of scattered light to the phase curve, the offset must become more pronounced to affect the phase curve and produce a visible asymmetry. Furthermore, for low albedos (e.g., high temperatures and low scattering contribution), inefficient heat recirculation is required to produce a phase curve at all. Upon increasing the scattering albedo, higher values of ϵ are allowed, but that reaches a maximum. Beyond this maximum, scattered light will dominate the phase curve, and again, ϵ must decrease to produce a significant thermal asymmetry (i.e., large daynight temperature differences).
Figure D.5 (left panel) illustrates a degeneracy in the “standard + asy” scenario, between A_{S}, d_{S}, l_{start} and l_{end}. Since the asymmetric phase curve requires a lower posteclipse amplitude to fit the data, the “evening” side must be brighter than the “morning” side. This can be achieved in two ways: either high A_{S} and correspondingly d_{S}< 1 and l_{start} and l_{end} delimiting part of the “morning” side, or low A_{S} and correspondingly d_{S}> 1 and l_{start} and l_{end} delimiting part of the evening side. However, note that from a physical standpoint, it is unclear how the albedo could be higher on the evening side. Mostly, it is assumed that clouds are responsible for the scattering. These are supposed to dissipate over time while circulating over the dayside hemisphere, therefore posteclipse maxima are not generally attributed to clouds (e.g., Demory et al. 2013; Esteves et al. 2015). Another possibility for posteclipse maxima being due to albedo changes would be the photodissociation of absorbers such as TiO or VO. In atomic form, the absorption would be much less efficient, hence the planetary albedo would increase. Investigating this possibility is, however, beyond the scope of this work.
Fig. 15 Joint credibility regions of recirculation and geometric albedo for TrES2b in the “standard” (black dots) and “standard + off” (blue dots) scenarios. Green contour: 1σ uncertainty region in Schwartz & Cowan (2015). Both this and previous work are consistent with each other. 
Figure 15 shows the constraints on recirculation and geometric albedo as inferred from our “standard” and “standard + off” scenarios, compared to results by Schwartz & Cowan (2015). As for CoRoT1b, we use joint credibility regions to illustrate our inferred range in the A_{G}ϵ plane. Both our results and the results of Schwartz & Cowan (2015) are consistent with each other and certainly agree better than for CoRoT1b (see Fig. 10). However, note that the 1σ region of Schwartz & Cowan (2015) contains roughly 30% of the points of the “standard” model, but only about 8% of the points of the “standard + off” model.
Similarly to CoRoT1b, we use the calculated A_{S} values to put constraints on the overall Bond albedo of TrES2b (Eqs. (36) and (37)). Using a_{V} ≈ 0.4, we obtain, based on the optical phase curve, A_{B}< 0.6. Putting together the dayside brightness temperature measurements from IRAC and Ks bands (O’Donovan et al. 2010; Croll et al. 2010), we then derive A_{B}< 0.68. These constraints are somewhat tighter than the ones derived for CoRoT1b because the spectral coverage is larger and TrES2b is a cooler planet.
5.4. HATP7b
Figure 16 shows our bestfit models of the different MCMC scenarios. In contrast to TrES2b, the phase curve of HATP7b is dominated by reflected light, rather than by the ellipsoidal variations (even though these are still clearly visible). Again, for clarity, the “free A” and “free T” scenarios are not shown. In Table D.3, we state bestfit parameters as well as 95% credibility regions for the parameters. Parameter posterior distributions are shown in Figs. D.6–D.8.
Fig. 16 HATP7b phase curve: comparison of bestfit models with data (red) and fit by Esteves et al. (2013; orange). Primary transit and secondary eclipse not shown. 
Figure 17 shows the marginalized posterior distributions for the inferred planetary mass. Also shown are results of previous phase curve modeling by Esteves et al. (2013) and Esteves et al. (2015) as well as RV mass determinations. Again, as for TrES2b, our estimated mass values are consistent with previous studies, and the determined planetary mass is not greatly affected by the choice of the phase curve model. Note that the formal uncertainties on planetary mass are somewhat smaller in our work than the RV uncertainties. This is mainly due to the excellent photometric quality of the phase curve and the fact that we fix stellar parameters, i.e., the stellar mass does not contribute to the final uncertainty on mass estimates. Note also the strong disagreement between mass estimates from Esteves et al. (2013) and Esteves et al. (2015; plain and dashed gray lines in Fig. 17, respectively). This is because the former uses separate beaming and ellipsoidal amplitudes as fitting parameters, while the latter uses planetary mass as a fitting parameter, as we do here. The beaming amplitude is adjusted to account for the asymmetry of the phase curve, thus planetary mass estimates from beaming and ellipsoidal amplitudes do not agree (4.2 compared to 1.6 M_{J}, see Table 5 in Esteves et al. 2013). Hence, it is clearly demonstrated that a separate fitting of both amplitudes can potentially lead to incorrect mass estimates.
Fig. 17 Marginalized posterior distributions for HATP7b M_{P} in different models. Mass from previous studies (including RV measurements) in gray. 
Figure 18 shows the inferred scattering albedo for the different scenarios. The values are broadly consistent with the values from Esteves et al. (2015) who find a geometric albedo of A_{G} ≈ 0.2, close to our values (recall , Eq. (8)). The fact that A_{S} is mostly independent of the specific planetary scenario suggests that the estimated value of A_{S} (0.26 <A_{S}< 0.34 at 95% confidence) is robust. The arithmetic mean for the combined scenarios is A_{S} = 0.28.
Fig. 18 Marginalized posterior distributions for HATP7b A_{S} in different models. A_{S} is mostly independent of the adopted model. 0.26 <A_{S}< 0.34 at 95% confidence. 
Fig. 19 Marginalized posterior distributions for HATP7b ϵ in different models. ϵ is strongly dependent on the adopted planetary scenario. 
As before, ϵ depends on the choice of the thermal model, as illustrated in Fig. 19. Similar to what has been found for TrES2b, the “standard + off” scenario requires a very different distribution for ϵ in order to produce a thermal offset in the phase curve. The fact that, for HATP7b, the distribution of the “standard + both” scenario is closer to the “standard + asy” distribution, suggests that for HATP7b, the asymmetry in the phase curve is better explained by scattered light than by thermal emission. This is supported by the BIC values (see Table D.3). Both the “standard +asy” and the “standard + both” scenarios are strongly favored by a probability of about 10^{5} compared to the “standard + off” scenario (ΔBIC ≈ 20). This result indicates that the preferred model explanation for the asymmetry is reflected light, rather than a thermal offset (but see above for a discussion of the physical problems of this solution). Our results quite clearly suggest asymmetric models rather than symmetric ones (ΔBIC > 400), again confirming previous phase curve analysis (Esteves et al. 2015).
The values in Table D.3 are somewhat high (3.7 for the preferred model). Usually, such a high value might indicate that either the model is not capturing correctly the physical behavior of the system, or that the errors are underestimated. However, the is dominated by a few points (7 outliers contribute 50% of the total ), and one particular point contributes around 15%. Since this is found for all fit scenarios (i.e., always the same outliers), this suggests that the error bars are indeed underestimated. When removing the apparently systematic outliers, the is reduced to less than 2. This in turn suggests a relatively good fit.
Figure D.7 (right panel) shows the correlations for the “standard + off” scenario, as discussed above for TrES2b. These correlations are much stronger and cleaner in this case, since the signaltonoise ratio of the phase curve is much better for HATP7b.
Figure 20 shows the constraints on recirculation and geometric albedo as inferred from our scenarios, compared to Schwartz & Cowan (2015). As above, we use joint credibility regions to illustrate our inferred range in the A_{G}ϵ plane. Note that the “standard+asy” model results and the lowalbedo part of the “standard+both” model overlap (see also Fig. D.8). It is clearly seen that our results and the results of Schwartz & Cowan (2015) strongly disagree, as was the case for CoRoT1b. However, because of the good data quality of the HATP7b phase curve, the disagreement is stronger than for CoRoT1b. There is no overlap between our 90% joint credibility regions and the 1σ uncertainty region of Schwartz & Cowan (2015). This is because of the very wellconstrained scattering albedo, which is nearly independent of the thermal and asymmetry model that was chosen (see also Fig. 18).
Fig. 20 Joint credibility regions of recirculation and geometric albedo for HATP7b in different scenarios. Green contour: 1σ uncertainty region in Schwartz & Cowan (2015). Both our and previous work strongly disagree on the inferred ϵ and A_{G} values. 
We use the calculated A_{S} values and measured IR dayside emission spectrum to constrain A_{B} for HATP7b (Eqs. (36) and (37)). As for TrES2b, we have a_{V} ≈ 0.4. Therefore, based on the optical phase curve, we find 0.11 <A_{B}< 0.72. The observed Spitzer spectrum translates into A_{B}< 0.87. These constraints are very loose because HATP7b is a rather hot planet, compared to TrES2b. NearIR measurements covering the 1–3 μm range (close to the Wien peak of the thermal radiation) would be highly desirable to further constrain the energy balance and add more constraints on the Bond albedo.
6. Discussion
For solar system objects, phase curves have been incredibly useful to determine the scattering properties of atmospheric particles (size distribution, vertical location and extent, composition). Groundbased data as well as spacecraft observations (e.g., Venus Express, Voyager 1 and 2, Pioneer 10 and 11) have been used to investigate, e.g., Venus (e.g., Arking & Potter 1968; García Muñoz et al. 2014; Petrova et al. 2015), Titan (e.g., Rages et al. 1983), Jupiter (e.g., Tomasko et al. 1978; Smith & Tomasko 1984), Saturn (e.g., Tomasko & Doose 1984) or Uranus (Rages et al. 1991; Pryor et al. 1997).
Large differences are found in the broadband phase curves of, e.g., Jupiter and Saturn (e.g., Dyudina et al. 2005) or Mars, Mercury and Venus (e.g., Mallama 2009). These are of course attributable to differences in cloud structure and composition, the absence or presence of an atmosphere, topographic surface features or the amount of dustcovered or bare regolith, to name but a few factors influencing the phase curves.
Sophisticated radiative transfer models in combination with cloud and aerosol models are needed to interpret these observations correctly and retrieve scattering properties.
In comparison, exoplanet studies suffer from the incredibly crude data available at present (in terms of signaltonoise ratio, spectral resolution or spectral coverage). Even though recent progress has been astonishing, we do not expect exoplanet data to approach SolarSystem quality in the near future. Hence, interpretation of exoplanet observations does not require models of comparable complexity yet, although more complex models, which take into account, for example, cloud formation or temperature gradients, have recently been published (e.g., Webber et al. 2015; Hu et al. 2015) and applied to, e.g., the wellcharacterized phase curve of Kepler7b. However, most studies, including this work, rely on simpler models and make strong assumptions to infer planetary and atmospheric properties.
For example, the model used in this work relies on the following two assumptions, in line with previous studies (e.g., Snellen et al. 2009; Cowan & Agol 2011; Schwartz & Cowan 2015):

Day and night hemispheres are assumed to be respectivelydescribed by a single, uniform temperature, without any longitudinal or latitudinal gradients. In reality, this is unlikelyto be true. Secondaryeclipse mapping of the hot JupiterHD 189733b has already demonstratedthat the brightness distribution is far from uniform (e.g., de Witet al. 2012). This can be interpreted as a nonuniform temperature distribution. IR phasecurves also clearly show temperature gradients (e.g., Knutsonet al. 2007, 2009;Crossfield et al. 2010). In the hypothetical norecirculation limit (ϵ = 0), nonuniformity effects mightproduce a thermal beaming dominated by the substellar pointthat could potentially enhance planetary emitted radiation (e.g.,Selsis et al. 2011; Schwartz &Cowan 2015). However, given the relatively lowsignaltonoise ratios, the number of effects contributing to theoptical phase curve and the large bandpass of Kepler and CoRoT,the optical phase curve is not expected to be very sensitive to thetemperature distribution.

The observed brightness temperature in a given spectral bandpass equals the bolometric equilibrium temperature hence constrains the energy budget of the atmosphere and can be related to the Bond albedo. This is a fundamental assumption that is unlikely to hold once better spectral resolution becomes available. As suggested by, e.g., Barclay et al. (2012), the photospheres for optical and IR observations (as well as for day and nightside emission) are probably located at different pressures. Hence, the observations would probe different temperatures and dynamical regimes. Depending on pressure, circulation and temperature regimes can be quite different (e.g., Parmentier et al. 2013; Agúndez et al. 2014; Showman et al. 2015). Hence, observed brightness temperatures in either spectral domain would not be necessarily related to the bolometric equilibrium temperature.
It is possible that these assumptions are not violated (or at least, not strongly) for many exoplanets and that they more or less hold. Our results, however, imply that optical and IR data lead to different conclusions for the same objects (in two out of three cases) when applying these assumptions. Therefore, it seems that they are too strong and overly simplified. It is a subject of future research to reconcile this finding with the current data quality, which does not necessarily warrant complex models or a level of sophistication much higher than the models presented here or in previous work.
We point out, however, that in the case for CoRoT1b, both our model results and the results by Schwartz & Cowan (2015) are marginally compatible, since their 1σ uncertainty regions and our 90% credibility regions slightly overlap. Therefore, an analysis of the CoRoT1b phase curve with the newly released, improved data pipeline might reduce the photometric uncertainties and provide a more decisive answer to resolve the apparent contradiction between phase curve and secondary eclipse analyses.
7. Conclusions
We have presented a simple, yet physically consistent, model of optical phase curves for exoplanets. It includes Lambertian scattering, thermal emission (under the assumption of uniform hemispheric temperatures), ellipsoidal variations and Doppler boosting. It can account for asymmetric phase curves by longitudinally asymmetric scattering albedos and an offset of thermal radiation compared to the substellar point.
This model has been used to analyze published phasecurve data of CoRoT1b, TrES2b and HATP7b. Results are then compared to an analysis of secondaryeclipse data of these planets by Schwartz & Cowan (2015).
We have shown that for CoRoT1b and HATP7b, inferred albedo and heat recirculation values from optical phase curves are different compared to previously published results. For TrES2b, both methods yield similar results.
We find that CoRoT1b has a rather higher scattering albedo than previously found. We find 0.11 <A_{S}< 0.3 at 95% confidence, which is in slight contrast with previous analyses, which found A_{S}< 0.15 (Snellen et al. 2009; Schwartz & Cowan 2015). Also, full phase curve analysis favors a strong redistribution of stellar incident energy to the nightside, contrary to previous studies, which suggested a very inefficient recirculation (Snellen et al. 2009; Schwartz & Cowan 2015). These contradictions are mainly because previous optical phase curve analysis of CoRoT1b by Snellen et al. (2009) considered only thermal emission.
In line with previous studies on the optical phase curve of HATP7b (e.g., Esteves et al. 2013, 2015), we find an appreciable albedo (A_{S} ≈ 0.3), slightly higher than inferred from secondary eclipse data. In contrast to previous studies based on secondary eclipse data (Schwartz & Cowan 2015), the analysis of the optical phase curve favors moderate to efficient heat recirculation. Asymmetric models are found to best fit the observed phase curve.
These differences between secondary eclipse and optical phase curve analyses occur most likely because optical and IR observations probe different atmospheric layers. Furthermore, our results suggest that some of the assumptions made (specifically, that observed brightness temperatures constrain the energy budget) are probably too strong and should be relaxed.
Future work will aim, among others, at reanalyzing further planets with published optical phase curves and reconciling different observations in the optical and the IR for CoRoT1b.
Note that we calculated the star’s solid angle, as seen from the planet, in Eq. (10) as . This assumes that R_{∗} ≪ r, a condition that is on the verge of breaking down for closein planets such as the ones considered in this work. However, detailed modeling (not shown here), which took the large angular extent of the star (up to tens of degrees) into account, showed little to no influence on resulting phase curves. Therefore, we retain our simplifying assumption in the following.
Retrieved from http://keplergo.arc.nasa.gov/ Instrumentation.shtml
Note that for tidallylocked planets, such as assumed here, the star does not move across the celestial sphere for an observer on the planet. Therefore, “morning” (i.e., sunrise) and “evening” (i.e., sunset) as such don’t exist. Rather, for eastward circulation, “morning” is defined as the terminator over which an air parcel would enter the dayside from the nightside. For illustration purposes, we will retain “morning” and “evening” throughout the text.
Acknowledgments
This study has received financial support from the French State in the frame of the “Investments for the future” Programme IdEx Bordeaux, reference ANR10IDEX0302. P.G. acknowledges support from the ERC Starting Grant (3DICE, grant agreement 336474). We thank the anonymous referee for a very constructive and positive feedback. Computer time for this study was provided in parts by the computing facilities MCIA (Mésocentre de Calcul Intensif Aquitain) of the Université de Bordeaux and of the Université de Pau et des Pays de l’Adour.
References
 Agúndez, M., Venot, O., Iro, N., et al. 2012, A&A, 548, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Agúndez, M., Parmentier, V., Venot, O., Hersant, F., & Selsis, F. 2014, A&A, 564, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Alonso, R., Alapini, A., Aigrain, S., et al. 2009, A&A, 506, 353 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Arking, A., & Potter, J. 1968, J. Atmosph. Sci., 25, 617 [NASA ADS] [CrossRef] [Google Scholar]
 Arras, P., & Socrates, A. 2010, ApJ, 714, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Barclay, T., Huber, D., Rowe, J. F., et al. 2012, ApJ, 761, 53 [NASA ADS] [CrossRef] [Google Scholar]
 Barge, P., Baglin, A., Auvergne, M., et al. 2008, A&A, 482, L17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Borucki, W. J., Koch, D., Jenkins, J., et al. 2009, Science, 325, 709 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Castelli, F., & Kurucz, R. L. 2004, ArXiv eprints [arXiv:astroph/0405087] [Google Scholar]
 Christiansen, J. L., Ballard, S., Charbonneau, D., et al. 2010, ApJ, 710, 97 [NASA ADS] [CrossRef] [Google Scholar]
 Claret, A., & Bloemen, S. 2011, A&A, 529, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Collier Cameron, A., Horne, K., Penny, A., & Leigh, C. 2002, MNRAS, 330, 187 [NASA ADS] [CrossRef] [Google Scholar]
 Cowan, N. B., & Agol, E. 2011, ApJ, 729, 54 [NASA ADS] [CrossRef] [Google Scholar]
 Croll, B., Albert, L., Lafreniere, D., Jayawardhana, R., & Fortney, J. J. 2010, ApJ, 717, 1084 [NASA ADS] [CrossRef] [Google Scholar]
 Crossfield, I. J. M., Hansen, B. M. S., Harrington, J., et al. 2010, ApJ, 723, 1436 [NASA ADS] [CrossRef] [Google Scholar]
 de Wit, J., Gillon, M., Demory, B.O., & Seager, S. 2012, A&A, 548, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Deming, D., Knutson, H., Agol, E., et al. 2011, ApJ, 726, 95 [NASA ADS] [CrossRef] [Google Scholar]
 Demory, B.O., de Wit, J., Lewis, N., et al. 2013, ApJ, 776, L25 [NASA ADS] [CrossRef] [Google Scholar]
 Désert, J.M., VidalMadjar, A., Lecavelier Des Etangs, A., et al. 2008, A&A, 492, 585 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dyudina, U. A., Sackett, P. D., Bayliss, D. D. R., et al. 2005, ApJ, 618, 973 [NASA ADS] [CrossRef] [Google Scholar]
 Esteves, L. J., De Mooij, E. J. W., & Jayawardhana, R. 2013, ApJ, 772, 51 [NASA ADS] [CrossRef] [Google Scholar]
 Esteves, L. J., De Mooij, E. J. W., & Jayawardhana, R. 2015, ApJ, 804, 150 [NASA ADS] [CrossRef] [Google Scholar]
 Faigler, S., & Mazeh, T. 2015, ApJ, 800, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Faigler, S., TalOr, L., Mazeh, T., Latham, D. W., & Buchhave, L. A. 2013, ApJ, 771, 26 [NASA ADS] [CrossRef] [Google Scholar]
 Feigelson, E., & Babu, G. 2012, Modern Statistical Methods for Astronomy (Cambridge University Press) [Google Scholar]
 ForemanMackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [CrossRef] [Google Scholar]
 García Muñoz, A., PérezHoyos, S., & SánchezLavega, A. 2014, A&A, 566, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gillon, M., Demory, B., Triaud, A. H. M. J., et al. 2009, A&A, 506, 359 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Goggans, P. M., & Chi, Y. 2004, AIP Conf. Ser., 707, 59 [NASA ADS] [CrossRef] [Google Scholar]
 Goodman, J., & Weare, J. 2010, Commun. Appl. Math. Comput. Sci., 5 [Google Scholar]
 Heng, K., & Showman, A. 2015, Ann. Rev. Earth Planet. Sciences, 43 [Google Scholar]
 Hou, F., Goodman, J., Hogg, D. W., Weare, J., & Schwab, C. 2012, ApJ, 745, 198 [NASA ADS] [CrossRef] [Google Scholar]
 Hu, R., Demory, B.O., Seager, S., Lewis, N., & Showman, A. P. 2015, ApJ, 802, 51 [NASA ADS] [CrossRef] [Google Scholar]
 Kane, S. R., & Gelino, D. M. 2010, ApJ, 724, 818 [NASA ADS] [CrossRef] [Google Scholar]
 Kipping, D. M., & Spiegel, D. S. 2011, MNRAS, 417, L88 [NASA ADS] [Google Scholar]
 Knutson, H. A., Charbonneau, D., Allen, L. E., et al. 2007, Nature, 447, 183 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Knutson, H. A., Charbonneau, D., Cowan, N. B., et al. 2009, ApJ, 690, 822 [NASA ADS] [CrossRef] [Google Scholar]
 Lecavelier Des Etangs, A., Pont, F., VidalMadjar, A., & Sing, D. 2008, A&A, 481, L83 [NASA ADS] [CrossRef] [Google Scholar]
 Loeb, A., & Gaudi, B. S. 2003, ApJ, 588, L117 [NASA ADS] [CrossRef] [Google Scholar]
 Madhusudhan, N., & Burrows, A. 2012, ApJ, 747, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Madhusudhan, N., Knutson, H., Fortney, J. J., & Barman, T. 2014, Protostars and Planets VI, 739 [Google Scholar]
 Mallama, A. 2009, Icarus, 204, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Mazeh, T., & Faigler, S. 2010, A&A, 521, L59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Meeus, J. 1991, Astronomical Algorithms, 2nd edn. (WillmannBell, Inc.) [Google Scholar]
 Mislis, D., Heller, R., Schmitt, J. H. M. M., & Hodgkin, S. 2012, A&A, 538, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 O’Donovan, F. T., Charbonneau, D., Mandushev, G., et al. 2006, ApJ, 651, L61 [NASA ADS] [CrossRef] [Google Scholar]
 O’Donovan, F. T., Charbonneau, D., Harrington, J., et al. 2010, ApJ, 710, 1551 [NASA ADS] [CrossRef] [Google Scholar]
 Pál, A., Bakos, G. Á., Torres, G., et al. 2008, ApJ, 680, 1450 [NASA ADS] [CrossRef] [Google Scholar]
 Parmentier, V., Showman, A. P., & Lian, Y. 2013, A&A, 558, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Penndorf, R. 1957, J. Opt. Soc. America (19171983), 47, 176 [Google Scholar]
 PerezBecker, D., & Showman, A. P. 2013, ApJ, 776, 134 [NASA ADS] [CrossRef] [Google Scholar]
 Petrova, E. V., Shalygina, O. S., & Markiewicz, W. J. 2015, Planet. Space Science, 113, 120 [NASA ADS] [CrossRef] [Google Scholar]
 Pfahl, E., Arras, P., & Paxton, B. 2008, ApJ, 679, 783 [NASA ADS] [CrossRef] [Google Scholar]
 Placek, B., Knuth, K. H., & Angerhausen, D. 2014, ApJ, 795, 112 [NASA ADS] [CrossRef] [Google Scholar]
 Pryor, W. R., West, R. A., & Simmons, K. E. 1997, Icarus, 127, 508 [NASA ADS] [CrossRef] [Google Scholar]
 Quintana, E. V., Rowe, J. F., Barclay, T., et al. 2013, ApJ, 767, 137 [NASA ADS] [CrossRef] [Google Scholar]
 Rages, K., Pollack, J. B., & Smith, P. H. 1983, J. Geophys. Res., 88, 8721 [NASA ADS] [CrossRef] [Google Scholar]
 Rages, K., Pollack, J. B., Tomasko, M. G., & Doose, L. R. 1991, Icarus, 89, 359 [NASA ADS] [CrossRef] [Google Scholar]
 Rauscher, E., & Kempton, E. M. R. 2014, ApJ, 790, 79 [NASA ADS] [CrossRef] [Google Scholar]
 Rogers, J. C., Apai, D., LopezMorales, M., Sing, D. K., & Burrows, A. 2009, ApJ, 707, 1707 [NASA ADS] [CrossRef] [Google Scholar]
 Rowe, J. F., Matthews, J. M., Seager, S., et al. 2006, ApJ, 646, 1241 [NASA ADS] [CrossRef] [Google Scholar]
 Rowe, J. F., Matthews, J. M., Seager, S., et al. 2008, ApJ, 689, 1345 [NASA ADS] [CrossRef] [Google Scholar]
 Schlawin, E., Zhao, M., Teske, J. K., & Herter, T. 2014, ApJ, 783, 5 [NASA ADS] [CrossRef] [Google Scholar]
 Schwartz, J. C., & Cowan, N. B. 2015, MNRAS, 449, 4192 [NASA ADS] [CrossRef] [Google Scholar]
 Selsis, F., Wordsworth, R. D., & Forget, F. 2011, A&A, 532, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Shardanand & Rao, A. D. P. 1977, NASA Technical Note D8442 [Google Scholar]
 Showman, A. P., & Guillot, T. 2002, A&A, 385, 166 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Showman, A. P., & Polvani, L. M. 2011, ApJ, 738, 71 [NASA ADS] [CrossRef] [Google Scholar]
 Showman, A. P., Lewis, N. K., & Fortney, J. J. 2015, ApJ, 801, 95 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, P. H., & Tomasko, M. G. 1984, Icarus, 58, 35 [NASA ADS] [CrossRef] [Google Scholar]
 Snellen, I. A. G., de Mooij, E. J. W., & Albrecht, S. 2009, Nature, 459, 543 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Spiegel, D. S., & Burrows, A. 2010, ApJ, 722, 871 [NASA ADS] [CrossRef] [Google Scholar]
 Stevenson, K. B., Désert, J.M., Line, M. R., et al. 2014, Science, 346, 838 [NASA ADS] [CrossRef] [Google Scholar]
 Sudarsky, D., Burrows, A., & Pinto, P. 2000, ApJ, 538, 885 [NASA ADS] [CrossRef] [Google Scholar]
 Tomasko, M. G., & Doose, L. R. 1984, Icarus, 58, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Tomasko, M. G., West, R. A., & Castillo, N. D. 1978, Icarus, 33, 558 [NASA ADS] [CrossRef] [Google Scholar]
 Van Eylen, V., Kjeldsen, H., ChristensenDalsgaard, J., & Aerts, C. 2012, Astron. Nachr., 333, 1088 [NASA ADS] [CrossRef] [Google Scholar]
 Van Eylen, V., Lindholm Nielsen, M., Hinrup, B., Tingley, B., & Kjeldsen, H. 2013, ApJ, 774, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Webber, M. W., Lewis, N. K., Marley, M., et al. 2015, ApJ, 804, 94 [NASA ADS] [CrossRef] [Google Scholar]
 Welsh, W. F., Orosz, J. A., Seager, S., et al. 2010, ApJ, 713, L145 [NASA ADS] [CrossRef] [Google Scholar]
 Williams, D. M., & Pollard, D. 2002, Int. J. Astrobiol., 1, 61 [NASA ADS] [CrossRef] [Google Scholar]
 Winn, J. N., Johnson, J. A., Albrecht, S., et al. 2009, ApJ, 703, L99 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Model verification
We verify that the implementation of model equations leads to the correct limits for reflected and emitted light.
Equations (A.1) and (A.2) show the phase function for a standard Lambertian sphere. These equations have been used in most studies of optical phase curves so far.
Note that α_{0} is 0 at opposition and π at primary transit, contrary to α (see Eq. (2)).
Fig. A.1 Model test: reflected flux compared to exact Lambertian sphere. 
For thermal radiation, the phase function is described by the illuminated fraction L of the planetary disk, i.e., the dayside: (A.3)In Fig. A.1, we show the difference between the exact formulation (Eqs. (A.1) and (A.2)) and our model for the reflected component. The considered case is a hot Jupiter in a 10day orbit around a Sunlike star, at varying orbital inclinations. The amplitude of the signal is of the order of a few ppm (10^{6}). The difference is about three orders of magnitude less, which clearly indicates that the model correctly incorporates Lambertian scattering. The highfrequency structure in the residuals is due to the spatial discretization of the numerical model and has no effect on physical results.
Fig. A.2 Model test: Emitted flux compared to exact solution. 
Figure A.2 shows the comparison of our model to the exact solution for thermal radiation. In this case, we show a hot Jupiter in a 2day orbit around a Sunlike star (A_{B} = 0.15, ϵ = 0), in order to get an appreciable signal. The difference at peak amplitude is about 0.01 ppm for a total amplitude of ≈1.6 ppm. This amounts to an error of less than 1%, which we deem acceptable. Again, the highfrequency structure in the residuals is due to the spatial discretization of the numerical model and the time resolution.
Appendix B: Phase curves of transiting planets
If the data is of high enough quality, in terms of signaltonoise ratio or time resolution, more and more parameters can be added to the fit. However, when analyzing phase curves of transiting planets, some parameters can be related to one another selfconsistently, by the shape of the primary transit. For instance, the transit depth of the primary transit directly yields the radius ratio k_{r} between planet and star: (B.1)Furthermore, for circular orbits, transit duration and transit shape can be related to the orbital inclination i and the projected starplanet separation k_{p} in units of stellar radii: (B.2)Assuming M_{P} ≪ M_{∗}, we can write Kepler’s 3rd law as follows: (B.3)Since the orbital period P_{orb} of transiting planets is usually known to within a few minutes or better, and k_{p} and k_{r} are mostly determined to an accuracy of better than 1%, it is possible to calculate the stellar mass and the planetary radius, given a stellar radius. Equation (B.3) yields, in this case, an analytic relation between stellar radius and stellar mass.
Such a relation is shown in Fig. B.1 (blue line) for HATP7, using a period of P_{orb} = 2.204 days (Pál et al. 2008) and k_{p} = 4.1512 (Esteves et al. 2013). Also shown are stellar parameters taken from Pál et al. (2008) who used highresolution spectroscopy and Van Eylen et al. (2012) who used asteroseismology (see Tables B.1 and B.2 for a compilation).
Stellar parameters for HATP7.
Planetary parameters for HATP7b.
Fig. B.1 Consistency between reported radius and mass determinations for the star HATP7, by using the determined orbital period P_{orb} = 2.204 days, a/R_{∗} values from Table B.2, blue line is Eq. (B.3). 
It is clearly seen that fixing stellar parameters at R_{∗} = 1.84 R_{S} and M_{∗} = 1.47 M_{S}, as done by Esteves et al. (2013), results in inconsistent system parameters. These then introduce a significant error in estimating the planetary mass from the ellipsoidal variations.
To illustrate the effects on planetary mass estimates, we performed inverse modeling of the HATP7b phase curve, adopting the “standard” scenario from Table 2, i.e., fitting for mass, albedo (A_{B} = A_{S}) and heat recirculation. Consistent models use the stellar parameters from Van Eylen et al. (2012), i.e., R_{∗} = 1.90 R_{S} and M_{∗} = 1.36 M_{S}, whereas the inconsistent models use R_{∗} = 1.84 R_{S} and M_{∗} = 1.47 M_{S}, as done in Esteves et al. (2013).
Figure B.2 shows the residuals ΔC = C_{con}−C_{incon} of the bestfit models. As is clearly seen, both scenarios result in virtually identical fits, which only differ by about 0.2 ppm (compared to the roughly 75 ppm amplitude of the phase curve). Furthermore, both bestfit models result in similar values of 10.8 and 10.97, respectively.
Fig. B.2 Residuals between consistent (stellar mass and radius from Van Eylen et al. 2012) and inconsistent (stellar mass and radius from Pál et al. 2008) bestfit models of the HATP7b optical phase curve. 
All parameters except planetary mass are not affected by the choice of stellar parameters. Figure B.3 shows the marginalized posterior distributions for the planetary mass, for both sets of stellar parameters. It is clear that the estimated planetary mass varies by as much as 30%. In case of inconsistent stellar parameters, the planetary mass is severely overestimated, compared to RV results. Therefore, we chose the stellar parameters stated in Van Eylen et al. (2012) since these are consistent with parameters deduced from primary transit analysis.
Fig. B.3 Marginalized posterior distributions for HATP7b M_{P}, for different adopted stellar parameters in the standard model. 
Appendix C: Optional phase function choices
Appendix C.1: Empirical solar system phase functions
A few previous studies (e.g., Collier Cameron et al. 2002; Kane & Gelino 2010) used an empirically derived phase function instead of the Lambert phase function in Eq. (A.1). This phase function was obtained from a fit to optical observations of Venus and Jupiter.
When fitting the phase curve of HATP7b with the empirical phase function, we obtain a geometric albedo of about A_{G} ≈ 0.2, consistent with previous estimates using the Lambertian approximation. However, as shown in Fig. C.1, the estimated planetary mass is far larger. Even when changing from a uniform prior to a Gaussian prior based on RV measurements (1.78 ± 0.08, Esteves et al. 2015), the fitted mass is greatly overestimated. Hence, our results suggest that Eq. (C.2) as a particular choice of phase function is probably not correct for HATP7b. Possible reasons to explain this include, e.g., the higher temperature (much higher than both Venus and Jupiter, for which this particular phase function was derived) and a consequently much different atmospheric chemistry. Also, cloud properties could play a role, since HATP7b most likely has some form of silicate or iron clouds (see, e.g., comparison of Kepler7b and Jupiter in Webber et al. 2015). Even for Jupiter and Saturn, cloud properties are thought to be responsible for the difference in observed phase functions (e.g., Dyudina et al. 2005).
Such an impact of the choice of the phase function on mass estimates has also been discussed by Mislis et al. (2012).
Fig. C.1 Constraints on geometric albedo and planetary mass, as derived from the standard model using the empirical phase function of Eq. (C.2). 
Appendix C.2: Rayleigh scattering
To investigate the influence of Rayleigh scattering, we incorporated H_{2} and He Rayleigh scattering in the model. These two species are thought to form the major constituents of gasgiant atmospheres.
The Rayleigh scattering cross sections of H_{2} and He are calculated as (C.3)where σ_{ray,i} of species i is given in cm^{2} per molecule, λ in μm and λ_{0,i} is a reference wavelength where σ_{0,i} has been measured.
This approach is used in many approximative treatments of Rayleigh scattering (see, e.g., Lecavelier Des Etangs et al. 2008). For the values of λ_{0} and σ_{0,i} in Eq. (C.3), measurements from Shardanand & Rao (1977) were used in this work, as tabulated in Table C.1.
The optical depth in an atmospheric layer j due to Rayleigh scattering, τ_{ray,j}, is obtained with the following equation: (C.4)where σ_{k}, C_{k,j} are the Rayleigh cross section and the column density of species k respectively. We calculate the column density as (C.5)where c_{k,j} is the volume mixing ratio of species k in layer j, g_{P} planetary gravity, μ_{atm} the mean molecular weight of the atmosphere (≈2 for H_{2}dominated atmospheres) and P_{k} is the layer pressure. The atmospheric layers are approximately spaced evenly in log P from the “surface” pressure P_{S} to 10^{4} bar.
From there, the total optical depth τ_{ray} for use in Eq. (C.8) is obtained by summing the optical depths of each layer from the surface to the model lid: (C.6)In Fig. C.2, the used Rayleigh scattering cross sections of H_{2} are compared to measurements reported in the literature (Shardanand & Rao 1977) as well as different approximations used in various models (Table 2 of Penndorf 1957; Lecavelier Des Etangs et al. 2008).
For H_{2}, the agreement with measurements is very good, again to within the stated error bars of Shardanand & Rao (1977). Also, the agreement with the approximation of Lecavelier Des Etangs et al. (2008) is very good. The comparison with the parametrization of H_{2} Rayleigh scattering using Penndorf (1957) data is less good.
Fig. C.2 Comparison of H_{2} Rayleigh scattering cross sections. Relative deviations in % between model and data sources (as indicated). Vertical lines show measurement uncertainties. gray horizontal line indicates 0% deviation. 
At the “bottom” of the atmosphere, at a prescribed “surface” pressure P_{S}, we impose a Lambertian surface with scattering albedo A_{S}. Hence, Eq. (7) is modified, (C.7)where ω is the singlescattering albedo (set to unity in the allscattering, zeroabsorption approximation used here), T_{ray} is the transmission along the optical path calculated as (C.8)with τ_{ray} the (zenith) optical depth due to Rayleigh scattering (calculated at the midpoint of the spectral interval considered). The value of τ_{ray} will depend critically on the choice of P_{S} (see above, Eqs. (C.4) and (C.6)). Φ_{ray} is the phase function of Rayleigh scattering (C.9)with φ_{o} the angle between observer and incoming stellar light, i.e., cosφ_{o} = s·o.
Figure C.3 shows the effect of Rayleigh scattering on the phase curve. Together with the standard Lambertian scattering approximation of Eq. (7), we show phase curves with varying values of P_{S}. As expected, with increasing P_{S}, hence increasing contribution of Rayleigh scattering to the reflected light, the phase curve changes. For P_{S} = 1 bar, the atmosphere starts to become visible, and at P_{S} = 10 bar, dominates over the “surface” contribution.
Fig. C.3 Effect of Rayleigh scattering on the phase curve, for different values of P_{S}. 1 ppm error bar to the left. See text for discussion. 
Atmospheric modeling of hot Jupiters predicts the formation of clouds around the 10^{2} bar layer or even at lower pressures (e.g., Parmentier et al. 2013; Webber et al. 2015). This implies that the reflecting “surface” is at P_{S} < 10^{2} bar. Furthermore, optical absorption by, e.g., alkali metals or TiO/VO, greatly increases with pressure. Based on cross sections and solar abundances presented by Désert et al. (2008), Fig. C.4 shows the transmission due to TiO and VO absorption assuming P_{S} = 10^{4} bar. Figure C.4 suggests that not much radiation is expected to penetrate to levels where Rayleigh scattering becomes important.
Fig. C.4 Transmission due to TiO and VO, as a function of wavelength. See text for details. 
Hence, we would expect that Rayleigh scattering does not play a large role, and we neglect it in our phasecurve studies (equivalent to setting P_{S} = 0). Therefore, phase curves are calculated with the Lambert approximation.
Appendix D: MCMC results
Fig. D.1 CoRoT1b phase curve constraints: posterior projections for the standard model (left) and noscattering model (right). Dashed vertical lines represent marginalized 95% credibility regions. Smoothed 68% and 95% credibility regions in dark and light gray shade, respectively. Note that the results of both models are almost contrary to each other. 
Fig. D.2 CoRoT1b phase curve constraints: posterior projections for the freealbedo model (left) and freetemperature model (right). Dashed vertical lines represent marginalized 95% credibility regions. Smoothed 68% and 95% credibility regions in dark and light gray shade, respectively. 
Fig. D.3 TrES2b phase curve constraints: posterior projections for the “standard” model (left) and “free A” model (right). Dashed vertical lines represent marginalized 95% credibility regions. Smoothed 68% and 95% credibility regions in dark and light gray shade, respectively. 
Fig. D.4 TrES2b phase curve constraints: posterior projections for the “free T” model (left) and “standard + off” model (right). Dashed vertical lines represent marginalized 95% credibility regions. Smoothed 68% and 95% credibility regions in dark and light gray shade, respectively. 
Fig. D.5 TrES2b phase curve constraints: posterior projections for the “standard + asy” model (left) and “standard + both” model (right). Dashed vertical lines represent marginalized 95% credibility regions. Smoothed 68% and 95% credibility regions in dark and light gray shade, respectively. 
Fig. D.6 HATP7b phase curve constraints: posterior projections for the “standard” model (left) and “free A” model (right). Dashed vertical lines represent marginalized 95% credibility regions. Smoothed 68% and 95% credibility regions in dark and light gray shade, respectively. 
Fig. D.7 HATP7b phase curve constraints: posterior projections for the “free T” model (left) and “standard + off” model (right). Dashed vertical lines represent marginalized 95% credibility regions. Smoothed 68% and 95% credibility regions in dark and light gray shade, respectively. 
Fig. D.8 HATP7b phase curve constraints: posterior projections for the “standard + asy” model (left) and “standard + both” model (right). Dashed vertical lines represent marginalized 95% credibility regions. Smoothed 68% and 95% credibility regions in dark and light gray shade, respectively. 
Appendix E: MCMC convergence diagnostics
Fig. E.1 Trace plots of model parameters for the CoRoT1b “standard” (left) and “no scattering” (right) scenarios. Blue line traces the ensemble median, red lines correspond to the [0.16, 0.84] percentiles. 
Fig. E.2 Trace plots of model parameters for the CoRoT1b “free alb” (left) and “free temp” (right) scenarios. Blue line traces the ensemble median, red lines correspond to the [0.16, 0.84] percentiles. 
Fig. E.3 Trace plots of model parameters for the TrES2b “standard” (left) and “free alb” (right) scenarios. Blue line traces the ensemble median, red lines correspond to the [0.16, 0.84] percentiles. 
Fig. E.4 Trace plots of model parameters for the TrES2b “asymmetric” (left) and “free temp” (right) scenarios. Blue line traces the ensemble median, red lines correspond to the [0.16, 0.84] percentiles. 
Fig. E.5 Trace plots of model parameters for the TrES2b “offset” (left) and “both” (right) scenarios. Blue line traces the ensemble median, red lines correspond to the [0.16, 0.84] percentiles. 
Fig. E.6 Trace plots of model parameters for the HATP7b “standard” (left) and “free alb” (right) scenarios. Blue line traces the ensemble median, red lines correspond to the [0.16, 0.84] percentiles. 
Fig. E.7 Trace plots of model parameters for the HATP7b “asymmetric” (left) and “free temp” (right) scenarios. Blue line traces the ensemble median, red lines correspond to the [0.16, 0.84] percentiles. 
Fig. E.8 Trace plots of model parameters for the HATP7b “offset” (left) and “both” (right) scenarios. Blue line traces the ensemble median, red lines correspond to the [0.16, 0.84] percentiles. 
GelmanRubin statistics for the CoRoT1b scenarios.
GelmanRubin statistics for the TrES2b scenarios.
GelmanRubin statistics for the HATP7b scenarios.
All Tables
All Figures
Fig. 1 Basic sketch of the geometry of the model. 

In the text 
Fig. 2 Effect of thermal emission on the phase curve. See text for discussion. 

In the text 
Fig. 3 Illustration of the asymmetric models. Left: reflective dayside (red) with bright “morning” (green). Right: thermal offset modeled as a shifted dayside (red). See text for discussion. 

In the text 
Fig. 4 Trace plots of model parameters for the HATP7b “standard + asy” model (see Table 2). Blue line traces the ensemble median, red lines correspond to the [0.16, 0.84] percentiles (the dark red region in Fig. 5 below). 

In the text 
Fig. 5 Example of determination of uncertainty ranges via the cumulative probability distribution: Scattering albedo of CoRoT1b in the “standard” scenario (see text for further details). 68% credibility region in dark red, 95% credibility region in light red. Dashed lines indicate maximum a posteriori (MAP) value. 

In the text 
Fig. 6 CoRoT1b redchannel phase curve: Comparison of bestfit models with data (red) and fit by Snellen et al. (2009; yellow). Primary transit and secondary eclipse not shown. 

In the text 
Fig. 7 Marginalized posterior distributions for CoRoT1b scattering albedo in different models. A_{S} is approximately independent of the thermal model (0.11 <A_{S}< 0.3 at 95% confidence in the “free A” model). 

In the text 
Fig. 8 Marginalized posterior distributions for CoRoT1b ϵ (left) and A_{B} (right) in different models. Constraints on A_{B} are weak. ϵ is unconstrained and depends on the thermal model. 

In the text 
Fig. 9 Marginalized posterior distributions for CoRoT1b dayside and nightside temperatures in the “free T” model. 

In the text 
Fig. 10 Joint credibility regions of recirculation and geometric albedo (see Eq. (8)) for CoRoT1b in the “standard” (red dots) and “no scattering” (blue dots) scenarios. Orange contour: 1σ uncertainty region in Schwartz & Cowan (2015). Our “standard” model strongly disagrees with previous work. 

In the text 
Fig. 11 TrES2b phase curve: comparison of bestfit models with data (red) and fit by Barclay et al. (2012; orange). Primary transit and secondary eclipse not shown. 

In the text 
Fig. 12 Marginalized posterior distributions for TrES2b A_{S} in different models. TrES2b is a dark planet in all models (A_{S}< 0.03 at 95% confidence). 

In the text 
Fig. 13 Marginalized posterior distributions for TrES2b ϵ in different models. ϵ depends strongly on the chosen model. 

In the text 
Fig. 14 Marginalized posterior distributions for TrES2b M_{P} in different models. Mass from previous studies (including RV measurements) in gray. 

In the text 
Fig. 15 Joint credibility regions of recirculation and geometric albedo for TrES2b in the “standard” (black dots) and “standard + off” (blue dots) scenarios. Green contour: 1σ uncertainty region in Schwartz & Cowan (2015). Both this and previous work are consistent with each other. 

In the text 
Fig. 16 HATP7b phase curve: comparison of bestfit models with data (red) and fit by Esteves et al. (2013; orange). Primary transit and secondary eclipse not shown. 

In the text 
Fig. 17 Marginalized posterior distributions for HATP7b M_{P} in different models. Mass from previous studies (including RV measurements) in gray. 

In the text 
Fig. 18 Marginalized posterior distributions for HATP7b A_{S} in different models. A_{S} is mostly independent of the adopted model. 0.26 <A_{S}< 0.34 at 95% confidence. 

In the text 
Fig. 19 Marginalized posterior distributions for HATP7b ϵ in different models. ϵ is strongly dependent on the adopted planetary scenario. 

In the text 
Fig. 20 Joint credibility regions of recirculation and geometric albedo for HATP7b in different scenarios. Green contour: 1σ uncertainty region in Schwartz & Cowan (2015). Both our and previous work strongly disagree on the inferred ϵ and A_{G} values. 

In the text 
Fig. A.1 Model test: reflected flux compared to exact Lambertian sphere. 

In the text 
Fig. A.2 Model test: Emitted flux compared to exact solution. 

In the text 
Fig. B.1 Consistency between reported radius and mass determinations for the star HATP7, by using the determined orbital period P_{orb} = 2.204 days, a/R_{∗} values from Table B.2, blue line is Eq. (B.3). 

In the text 
Fig. B.2 Residuals between consistent (stellar mass and radius from Van Eylen et al. 2012) and inconsistent (stellar mass and radius from Pál et al. 2008) bestfit models of the HATP7b optical phase curve. 

In the text 
Fig. B.3 Marginalized posterior distributions for HATP7b M_{P}, for different adopted stellar parameters in the standard model. 

In the text 
Fig. C.1 Constraints on geometric albedo and planetary mass, as derived from the standard model using the empirical phase function of Eq. (C.2). 

In the text 
Fig. C.2 Comparison of H_{2} Rayleigh scattering cross sections. Relative deviations in % between model and data sources (as indicated). Vertical lines show measurement uncertainties. gray horizontal line indicates 0% deviation. 

In the text 
Fig. C.3 Effect of Rayleigh scattering on the phase curve, for different values of P_{S}. 1 ppm error bar to the left. See text for discussion. 

In the text 
Fig. C.4 Transmission due to TiO and VO, as a function of wavelength. See text for details. 

In the text 
Fig. D.1 CoRoT1b phase curve constraints: posterior projections for the standard model (left) and noscattering model (right). Dashed vertical lines represent marginalized 95% credibility regions. Smoothed 68% and 95% credibility regions in dark and light gray shade, respectively. Note that the results of both models are almost contrary to each other. 

In the text 
Fig. D.2 CoRoT1b phase curve constraints: posterior projections for the freealbedo model (left) and freetemperature model (right). Dashed vertical lines represent marginalized 95% credibility regions. Smoothed 68% and 95% credibility regions in dark and light gray shade, respectively. 

In the text 
Fig. D.3 TrES2b phase curve constraints: posterior projections for the “standard” model (left) and “free A” model (right). Dashed vertical lines represent marginalized 95% credibility regions. Smoothed 68% and 95% credibility regions in dark and light gray shade, respectively. 

In the text 
Fig. D.4 TrES2b phase curve constraints: posterior projections for the “free T” model (left) and “standard + off” model (right). Dashed vertical lines represent marginalized 95% credibility regions. Smoothed 68% and 95% credibility regions in dark and light gray shade, respectively. 

In the text 
Fig. D.5 TrES2b phase curve constraints: posterior projections for the “standard + asy” model (left) and “standard + both” model (right). Dashed vertical lines represent marginalized 95% credibility regions. Smoothed 68% and 95% credibility regions in dark and light gray shade, respectively. 

In the text 
Fig. D.6 HATP7b phase curve constraints: posterior projections for the “standard” model (left) and “free A” model (right). Dashed vertical lines represent marginalized 95% credibility regions. Smoothed 68% and 95% credibility regions in dark and light gray shade, respectively. 

In the text 
Fig. D.7 HATP7b phase curve constraints: posterior projections for the “free T” model (left) and “standard + off” model (right). Dashed vertical lines represent marginalized 95% credibility regions. Smoothed 68% and 95% credibility regions in dark and light gray shade, respectively. 

In the text 
Fig. D.8 HATP7b phase curve constraints: posterior projections for the “standard + asy” model (left) and “standard + both” model (right). Dashed vertical lines represent marginalized 95% credibility regions. Smoothed 68% and 95% credibility regions in dark and light gray shade, respectively. 

In the text 
Fig. E.1 Trace plots of model parameters for the CoRoT1b “standard” (left) and “no scattering” (right) scenarios. Blue line traces the ensemble median, red lines correspond to the [0.16, 0.84] percentiles. 

In the text 
Fig. E.2 Trace plots of model parameters for the CoRoT1b “free alb” (left) and “free temp” (right) scenarios. Blue line traces the ensemble median, red lines correspond to the [0.16, 0.84] percentiles. 

In the text 
Fig. E.3 Trace plots of model parameters for the TrES2b “standard” (left) and “free alb” (right) scenarios. Blue line traces the ensemble median, red lines correspond to the [0.16, 0.84] percentiles. 

In the text 
Fig. E.4 Trace plots of model parameters for the TrES2b “asymmetric” (left) and “free temp” (right) scenarios. Blue line traces the ensemble median, red lines correspond to the [0.16, 0.84] percentiles. 

In the text 
Fig. E.5 Trace plots of model parameters for the TrES2b “offset” (left) and “both” (right) scenarios. Blue line traces the ensemble median, red lines correspond to the [0.16, 0.84] percentiles. 

In the text 
Fig. E.6 Trace plots of model parameters for the HATP7b “standard” (left) and “free alb” (right) scenarios. Blue line traces the ensemble median, red lines correspond to the [0.16, 0.84] percentiles. 

In the text 
Fig. E.7 Trace plots of model parameters for the HATP7b “asymmetric” (left) and “free temp” (right) scenarios. Blue line traces the ensemble median, red lines correspond to the [0.16, 0.84] percentiles. 

In the text 
Fig. E.8 Trace plots of model parameters for the HATP7b “offset” (left) and “both” (right) scenarios. Blue line traces the ensemble median, red lines correspond to the [0.16, 0.84] percentiles. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.