© ESO, 2016

1. Introduction

In recent years, exoplanetary science has developed from a detection-centered astronomical science towards a characterization-centered 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.

Radiative transfer and atmospheric modeling by, e.g., Sudarsky et al. (2000) predicts very low optical albedos for cloud-free 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 close-in 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 re-radiated 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 day-night 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; Perez-Becker & 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., Perez-Becker & 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 well-characterized IR and optical measurements, namely CoRoT-1b, TrES-2b and HAT-P-7b, 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 set-up 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, non-Lambertian phase functions and stellar parameter uncertainties.

2. Forward model

2.1. Geometry

The basic geometry set-up of the star-planet 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 close-in 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: $\mathit{\theta }\mathrm{=}{\mathrm{90}}^{\mathrm{◦}}\mathrm{-}\mathit{i,}$(1)i.e., an edge-on orbit has i = 90°, and a face-on 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): $\mathit{\phi }\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{\mathrm{360}}^{\mathrm{◦}}\mathrm{-}\mathit{\alpha }\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{,}$(2)with $\mathit{\alpha }\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\mathit{T}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{+}{\mathit{\omega }}_{\mathrm{P}}\mathit{,}$(3)where T is the true anomaly and ω_P the argument of periastron. The true anomaly is calculated from the eccentric anomaly E: $\tan \left(\frac{\mathit{T}}{\mathrm{2}}\right)\mathrm{=}\sqrt{\frac{\mathrm{1}\mathrm{+}\mathit{e}}{\mathrm{1}\mathrm{-}\mathit{e}}}\tan \left(\frac{\mathit{E}}{\mathrm{2}}\right)\mathit{,}$(4)with e eccentricity. E is determined by numerically solving Kepler’s equation: $\mathit{E}\mathrm{-}\mathit{e}\sin \mathrm{\left(}\mathit{E}\mathrm{\right)}\mathrm{=}\mathit{M,}$(5)with M mean anomaly. The solution is obtained with a publicly available Fortran procedure¹, based on a method described in Meeus (1991). The mean anomaly M is given by $\mathit{M}\mathrm{=}\mathrm{2}\mathit{\pi }\mathrm{·}\left(\mathit{x}\mathrm{-}\mathrm{\lfloor }\mathit{x}\mathrm{\rfloor }\right)\mathit{,}$(6)where $\mathit{x}\mathrm{=}\frac{\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{peri}}}{{\mathit{P}}_{\mathrm{orb}}}$, 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. 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 ${\mathit{F}}_{\mathrm{r}\mathit{,}\mathrm{o}\mathit{,}\mathrm{cell}}\mathrm{=}\cos {\mathit{z}}_{\mathrm{s}}\mathrm{·}{\mathit{F}}_{\mathrm{\ast }\mathit{,}\mathrm{p}}\mathrm{·}\frac{{\mathit{A}}_{\mathrm{S}}}{\mathit{\pi }}\mathrm{·}\cos {\mathit{z}}_{\mathrm{o}}\mathrm{·}\frac{\mathrm{\Delta }\mathit{S}}{{\mathit{d}}^{\mathrm{2}}}\mathit{,}$(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²), d the observer-planet distance and A_S the (potentially wavelength-dependent) planetary scattering albedo. A_S is assumed to be constant with time, hence neglecting any time-dependent 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: ${\mathit{A}}_{\mathrm{G}}\mathrm{=}\frac{\mathrm{2}}{\mathrm{3}}{\mathit{A}}_{\mathrm{S}}\mathit{.}$(8)The surface element ΔS of the cell, as seen from the planet’s center, is calculated as $\mathrm{\Delta }\mathit{S}\mathrm{=}\mathrm{\Delta \Omega }\mathrm{·}{\mathit{R}}_{\mathrm{p}}^{\mathrm{2}}\mathit{,}$(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: ${\mathit{F}}_{\mathrm{\ast }\mathit{,}\mathrm{p}}\mathrm{=}\mathit{\pi }{\left(\frac{{\mathit{R}}_{\mathrm{\ast }}}{\mathit{r}}\right)}^{\mathrm{2}}{\mathrm{\int }}_{{\mathit{\lambda }}_{\mathrm{low}}}^{{\mathit{\lambda }}_{\mathrm{high}}}{\mathit{I}}_{\mathrm{\ast }\mathit{,}\mathrm{s}}{\mathit{q}}_{\mathit{I}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{d}\mathit{\lambda ,}$(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 star-planet distance and R_∗ is the stellar radius³. 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)⁴. Instrumental filter functions T_l are taken, for CoRoT-1, from Snellen et al. (2009) and, for TrES-2 and HAT-P-7, the Kepler Handbook⁵.

The time-dependent planet-star distance r is calculated with $\mathit{r}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\mathit{a}\frac{\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}}{\mathrm{1}\mathrm{+}\mathit{e}\cos \mathrm{\left(}\mathit{T}\mathrm{\right(}\mathit{t}\mathrm{\left)}\mathrm{\right)}}\mathit{,}$(11)with a the semi-major axis.

The total reflected flux F_r,o received by the observer is the sum over all contributing cells: ${\mathit{F}}_{\mathrm{r}\mathit{,}\mathrm{o}}\mathrm{=}\sum _{\mathrm{(}{n}\mathrm{·}{s}\hspace{0.17em}⩾\hspace{0.17em}\mathrm{0}\mathrm{)}\mathrm{\wedge }\mathrm{(}{n}\mathrm{·}{o}\hspace{0.17em}⩾\hspace{0.17em}\mathrm{0}\mathrm{)}}{\mathit{F}}_{\mathrm{r}\mathit{,}\mathrm{o}\mathit{,}\mathrm{cell}}\mathit{.}$(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 long-period (and consequently colder) planets. However, for close-in 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 re-distribution ϵ. Temperatures are then calculated, based on Cowan & Agol (2011): ${\mathit{T}}_{\mathrm{day}}^{\mathrm{4}}\mathrm{=}{\mathit{T}}_{\mathrm{\ast }}^{\mathrm{4}}\frac{{\mathit{R}}_{\mathrm{\ast }}^{\mathrm{2}}}{{\mathit{a}}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{0.5}}}\mathrm{·}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{A}}_{\mathrm{B}}\mathrm{)}\left(\frac{\mathrm{2}}{\mathrm{3}}\mathrm{-}\frac{\mathrm{5}}{\mathrm{12}}\mathit{ϵ}\right)\mathrm{=}{\mathit{T}}_{\mathrm{0}}^{\mathrm{4}}\mathrm{·}\left(\frac{\mathrm{2}}{\mathrm{3}}\mathrm{-}\frac{\mathrm{5}}{\mathrm{12}}\mathit{ϵ}\right)\mathit{,}$(13)and ${\mathit{T}}_{\mathrm{night}}^{\mathrm{4}}\mathrm{=}{\mathit{T}}_{\mathrm{\ast }}^{\mathrm{4}}\frac{{\mathit{R}}_{\mathrm{\ast }}^{\mathrm{2}}}{{\mathit{a}}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{0.5}}}\mathrm{·}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{A}}_{\mathrm{B}}\mathrm{)}\frac{\mathit{ϵ}}{\mathrm{4}}\mathrm{=}{\mathit{T}}_{\mathrm{0}}^{\mathrm{4}}\mathrm{·}\frac{\mathit{ϵ}}{\mathrm{4}}\mathit{,}$(14)where T_∗ is the stellar effective temperature and the additional factor (1−e²)^0.5 accounts for the mean flux received over an orbit (Williams & Pollard 2002).

The parameter ϵ is a re-parameterization 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: ${\mathit{L}}_{\mathrm{P}}\mathrm{=}\frac{\mathrm{1}}{\mathit{f}}\mathrm{·}\mathit{\pi }{\mathit{R}}_{\mathrm{P}}^{\mathrm{2}}\mathit{\sigma }{\mathit{T}}_{\mathrm{day}}^{\mathrm{4}}\mathrm{=}\frac{\mathrm{1}}{\mathit{f}}{\mathit{F}}_{\mathrm{d}}\mathit{.}$(15)For perfect heat recirculation, the entire planet is at a uniform temperature, and thus $\mathit{f}\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}}$, 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), $\mathit{f}\mathrm{=}\frac{\mathrm{2}}{\mathrm{3}}$.

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 re-arrange Eq. (15) to yield $\mathit{f}\mathrm{=}\frac{{\mathit{F}}_{\mathrm{d}}}{{\mathit{L}}_{\mathrm{P}}}\mathrm{=}\frac{{\mathit{F}}_{\mathrm{d}}}{{\mathit{F}}_{\mathrm{S}}}\mathrm{=}\frac{{\mathit{T}}_{\mathrm{day}}^{\mathrm{4}}}{{\mathit{T}}_{\mathrm{0}}^{\mathrm{4}}}\mathrm{·}$(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 day-night 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; Perez-Becker & 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 ${\mathit{F}}_{\mathrm{t}\mathit{,}\mathrm{o}\mathit{,}\mathrm{cell}}\mathrm{\left(}{\mathit{T}}_{\mathrm{c}}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{{\mathit{\lambda }}_{\mathrm{low}}}^{{\mathit{\lambda }}_{\mathrm{high}}}\mathit{B}\mathrm{\left(}{\mathit{T}}_{\mathrm{c}}\mathit{,\lambda }\mathrm{\right)}{\mathit{q}}_{\mathrm{I}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{d}\mathit{\lambda }\mathrm{·}\cos {\mathit{z}}_{\mathrm{o}}\mathrm{·}\frac{\mathrm{\Delta }\mathit{S}}{{\mathit{d}}^{\mathrm{2}}}\mathit{,}$(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., ): ${\mathit{F}}_{\mathrm{t}\mathit{,}\mathrm{o}}\mathrm{=}\sum _{\mathrm{(}{n}\mathrm{·}{s}\hspace{0.17em}⩾\hspace{0.17em}\mathrm{0}\mathrm{)}\mathrm{\wedge }\mathrm{(}{n}\mathrm{·}{o}\hspace{0.17em}⩾\hspace{0.17em}\mathrm{0}\mathrm{)}}{\mathit{F}}_{\mathrm{t}\mathit{,}\mathrm{o}\mathit{,}\mathrm{cell}}\mathrm{\left(}{\mathit{T}}_{\mathrm{day}}\mathrm{\right)}\mathrm{+}\sum _{\mathrm{(}{n}\mathrm{·}{s}\hspace{0.17em}⩽\hspace{0.17em}\mathrm{0}\mathrm{)}\mathrm{\wedge }\mathrm{(}{n}\mathrm{·}{o}\hspace{0.17em}⩾\hspace{0.17em}\mathrm{0}\mathrm{)}}{\mathit{F}}_{\mathrm{t}\mathit{,}\mathrm{o}\mathit{,}\mathrm{cell}}\mathrm{\left(}{\mathit{T}}_{\mathrm{night}}\mathrm{\right)}\mathit{.}$(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.

Figure 2 shows the reflective and the thermal contribution for a hot Jupiter planet in a 2-day orbit around a Sun-like 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 (dot-dashed 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, non-uniform 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 star-planet 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:

$\begin{array}{ccc}& & \\ & & \end{array}$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): ${\mathit{A}}_{\mathrm{ell}}\mathrm{=}{\mathit{\alpha }}_{\mathrm{ell}}\frac{{\mathit{M}}_{\mathrm{p}}}{{\mathit{M}}_{\mathrm{\ast }}}{\left(\frac{{\mathit{R}}_{\mathrm{\ast }}}{\mathit{a}}\right)}^{\mathrm{3}}{\sin }^{\mathrm{2}}\mathrm{\left(}\mathit{i}\mathrm{\right)}\mathit{,}$(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: ${\mathit{\alpha }}_{\mathrm{ell}}\mathrm{=}\mathrm{0.15}\mathrm{·}\frac{\mathrm{(}\mathrm{15}\mathrm{+}\mathit{u}\mathrm{)}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{g}\mathrm{)}}{\mathrm{3}\mathrm{-}\mathit{u}}\mathrm{·}$(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 pre-calculated 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_doppis given by (see also Eq. (5) in Quintana et al. 2013): ${\mathit{A}}_{\mathrm{dopp}}\mathrm{=}\mathrm{(}\mathrm{3}\mathrm{-}{\mathit{\alpha }}_{\mathrm{dopp}}\mathrm{)}\frac{\mathit{K}}{\mathit{c}}\mathit{,}$(23)where c is the speed of light, K the radial velocity semi-amplitude 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: ${\mathit{\alpha }}_{\mathrm{dopp}}\mathrm{=}\frac{{\mathit{e}}^{\mathit{x}}\mathrm{(}\mathrm{3}\mathrm{-}\mathit{x}\mathrm{)}\mathrm{-}\mathrm{3}}{{\mathit{e}}^{\mathit{x}}\mathrm{-}\mathrm{1}}\mathit{,}$(24)where $\mathit{x}\mathrm{=}\frac{\mathit{h}\frac{\mathit{c}}{{\mathit{\lambda }}_{\mathrm{obs}}}}{{\mathit{k}}_{\mathrm{B}}{\mathit{T}}_{\mathrm{\ast }}}$ (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 TrES-2b (their work, 3.71, to 3.87 with our approach) and HAT-P-7b (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) $\mathit{K}\mathrm{=}{\left(\frac{\mathrm{2}\mathit{\pi G}}{{\mathit{P}}_{\mathrm{orb}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}\frac{{\mathit{M}}_{\mathrm{p}}}{{\mathit{M}}_{\mathrm{\ast }}^{\mathrm{2}\mathit{/}\mathrm{3}}}\sin \mathrm{\left(}\mathit{i}\mathrm{\right)}\frac{\mathrm{1}\mathrm{+}\mathit{e}\cos {\mathit{\omega }}_{\mathrm{P}}}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{0.5}}}\mathrm{·}$(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”⁶ 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 re-radiated. 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 reflected-light component of the phase curve, we implement a simple dark-bright 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.

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 time-dependent contrast C between star and planet is then $\mathit{C}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\frac{{\mathit{F}}_{\mathrm{R}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{+}{\mathit{F}}_{\mathrm{E}}\mathrm{\left(}\mathit{t}\mathrm{\right)}}{{\mathit{F}}_{\mathrm{\ast }\mathit{,}\mathrm{o}}}\mathrm{+}{\mathit{C}}_{\mathrm{ell}}\mathrm{+}{\mathit{C}}_{\mathrm{dopp}}\mathit{.}$(26)The stellar flux at the observer, F_{∗ ,o}, is calculated analogous to F_{∗ ,p} (see Eq. (10)): ${\mathit{F}}_{\mathrm{\ast }\mathit{,}\mathrm{o}}\mathrm{=}\mathit{\pi }{\left(\frac{{\mathit{R}}_{\mathrm{\ast }}}{\mathit{d}}\right)}^{\mathrm{2}}{\mathrm{\int }}_{{\mathit{\lambda }}_{\mathrm{low}}}^{{\mathit{\lambda }}_{\mathrm{high}}}{\mathit{I}}_{\mathrm{\ast }\mathit{,}\mathrm{s}}{\mathit{q}}_{\mathrm{I}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{d}\mathit{\lambda }\mathit{.}$(27)

2.7. Model parameter summary

In total, our full physical model as described above contains up to 19 parameters, listed in Table 1.

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

We use the Bayesian formalism to calculate posterior probability values p(V_P | D) for the parameter vector V_P in the model, given a set D of observations. $\mathit{p}\mathrm{\left(}{\mathit{V}}_{\mathrm{P}}\mathrm{\right|}\mathit{D}\mathrm{)}\mathrm{\propto }\mathit{p}\mathrm{(}\mathit{D}\mathrm{\left|}{\mathit{V}}_{\mathrm{P}}\mathrm{\right)}\mathrm{·}\mathit{p}\mathrm{\left(}{\mathit{V}}_{\mathrm{P}}\mathrm{\right)}\mathit{.}$(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 Foreman-Mackey 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 auto-correlation lengths for each parameter). The first few auto-correlation lengths were considered as burn-in and discarded for the calculation of parameter uncertainties. Convergence was checked by inspecting visually the evolution of the mean of the entire ensemble and calculating the Gelman-Rubin 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.

Uncertainty ranges are calculated by marginalizing over the posterior distribution, thinned by 60 steps, covering roughly 1–3 auto-correlation 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] median-centered percentiles, respectively, of the cumulative probability distributions (CDF). If the parameter distribution were to be Gaussian, these credibility regions would correspond to the 1 and 2σ uncertainties, respectively. Figure 5 illustrates our method to determine credibility regions and best-fit parameters. Best-fit 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).

3.2. Goodness-of-fit criteria and model comparison

We use three standard quantities (e.g., Feigelson & Babu 2012) to evaluate the best-fitting models obtained from the MCMC model: χ², the reduced ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}$ and the Bayesian Information Criterion (BIC). These are evaluated for the MAP parameter set (not the median of the CDF).

χ², the sum of the weighted squared residuals, is defined as ${\mathit{\chi }}^{\mathrm{2}}\mathrm{=}\sum _{\mathit{j}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{1}}^{{\mathit{N}}_{\mathrm{D}}}{\left(\frac{{\mathit{M}}_{{\mathit{V}}_{\mathrm{P}}\mathit{,j}}\mathrm{-}{\mathit{D}}_{\mathit{j}}}{{\mathit{\sigma }}_{\mathit{j}}}\right)}^{\mathrm{2}}\mathit{,}$(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 ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}$ is generally calculated by ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}\mathrm{=}\frac{{\mathit{\chi }}^{\mathrm{2}}}{{\mathit{N}}_{\mathrm{D}}\mathrm{-}{\mathit{N}}_{\mathrm{P}}}\mathit{,}$(30)where N_P is the number of parameters. In order to be considered a good fit, ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}$ should be of the order of unity.

Note that by adding more parameters to the fitting model, the χ² will decrease. Therefore, we use another criterion, the BIC (valid when N_D ≫ N_P, which is the case for our calculations), defined as $\mathrm{BIC}\mathrm{=}{\mathit{\chi }}^{\mathrm{2}}\mathrm{+}{\mathit{N}}_{\mathrm{P}}\ln \mathrm{\left(}{\mathit{N}}_{\mathrm{D}}\mathrm{\right)}\mathit{.}$(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 ${\mathit{p}}_{\mathrm{M}}\mathrm{=}{\mathit{e}}^{\mathrm{-}\frac{\mathrm{\Delta }\mathrm{BIC}}{\mathrm{2}}}\mathit{,}$(32)where ΔBIC = BIC_M1−BIC_M0 is the difference between the model M₁ under consideration, and the model M₀ that minimizes the BIC.

The BIC is an approximation to the evidence (or marginal likelihood) E_model of a given model: ${\mathit{E}}_{\mathrm{model}}\mathrm{=}\mathrm{\int }\mathit{p}\mathrm{\left(}\mathit{D}\mathrm{\right|}{\mathit{V}}_{\mathrm{P}}\mathrm{)}\mathrm{·}\mathit{p}\mathrm{(}{\mathit{V}}_{\mathrm{P}}\mathrm{)}\mathrm{d}{\mathit{V}}_{\mathrm{P}}\mathit{.}$(33)The calculation of the integral in Eq. (33) is in most case computationally very expensive. In emcee, this integral is estimated using a so-called thermodynamic integration, based on an algorithm proposed by Goggans & Chi (2004). However, the calculation is very time-consuming (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 set-up

4.1. Planets

4.1.1. CoRoT-1b

CoRoT-1b 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 CoRoT-1b 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 CoRoT-1b 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 red-channel 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 secondary-eclipse 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. TrES-2b

TrES-2b 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. Ground-based 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 day-night 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 secondary-eclipse 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. HAT-P-7b

As TrES-2b, HAT-P-7b 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 3-year 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 radial-velocity 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 CoRoT-1b in terms of the normalized flux F_norm, normalized to the primary, instead of secondary, eclipse (i.e., without transit, the flux at zero-phase would be unity). Therefore, we calculate the forward model as: ${\mathit{F}}_{\mathrm{norm}}\mathrm{\left(}\mathit{\alpha }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}\mathrm{+}\mathit{C}\mathrm{\left(}\mathit{\alpha }\mathrm{\right)}}{\mathrm{1}\mathrm{+}\mathit{C}\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{·}$(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, non-physical parameter, a zero-point flux offset f₀. This parameter is related to the data reduction and not a priori linked to any physical characteristics of the star-planet system. The light curve L_C is then described with the following equation: ${\mathit{L}}_{\mathrm{C}}\mathrm{\left(}\mathit{\alpha }\mathrm{\right)}\mathrm{=}\mathit{C}\mathrm{\left(}\mathit{\alpha }\mathrm{\right)}\mathrm{+}{\mathit{f}}_{\mathrm{0}}\mathit{.}$(35)We will use this equation for our analysis of TrES-2b and HAT-P-7b. Hence, f₀ 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 CoRoT-1).

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: ${\mathit{a}}_{\mathit{V}}\mathrm{·}{\mathit{A}}_{\mathrm{S}}\mathrm{\le }{\mathit{A}}_{\mathrm{B}}\mathrm{\le }{\mathit{a}}_{\mathit{V}}\mathrm{·}{\mathit{A}}_{\mathrm{S}}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{a}}_{\mathit{V}}\mathrm{)}\mathit{.}$(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 phase-curve analysis by Snellen et al. (2009), we then use a fourth model approach for CoRoT-1b. 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 CoRoT-1b assume symmetric phase curves and do not fit for planetary mass because of the relatively low signal-to-noise ratio and reduced phase resolution of the binned phase curve.

For TrES-2b and HAT-P-7b, 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 re-emitted on the day and night hemispheres. $\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{A}}_{\mathrm{B}}\mathrm{)}\mathrm{·}{\mathit{F}}_{\mathrm{\ast }\mathit{,}\mathrm{p}}\mathrm{=}\mathrm{2}\left(\sum _{\mathit{j}}^{{\mathit{N}}_{\mathrm{day}}}{\mathit{F}}_{\mathit{j,}\mathrm{day}}\mathrm{+}\sum _{\mathit{j}}^{{\mathit{N}}_{\mathrm{day}}}{\mathit{F}}_{\mathit{j,}\mathrm{night}}\right)\mathit{,}$(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 Gelman-Rubin (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 inter-chain variance is large compared to the variance of individual chains (in our cases, sometimes orders of magnitude). For strongly non-linearly correlated parameters (e.g., in the “standard + off” and “standard + both” scenarios of HAT-P-7b, 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 TrES-2b (see Figs. E.4 and E.5, Table E.2). Since the parameters describing the albedo asymmetry are tightly anti-correlated, 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. CoRoT-1b

Figure 6 shows the data and the best-fit model from Snellen et al. (2009) as well as our best-fit models from the various planetary scenarios (Table 2). Since the “free T” and the “free albedo” best-fit models are virtually identical, only the latter is shown, for clarity. Clearly, the standard and free-albedo best-fit 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.

Best-fit parameters as well as goodness-of-fit 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.

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 CoRoT-1b (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).

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 secondary-eclipse 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.

Our results are somewhat contrary to the results from the IR secondary-eclipse 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 CoRoT-1b is probably a low-albedo 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 two-dimensional 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 no-scattering case is equivalent to imposing a strong prior on the scattering albedo (A_S = 0) that seems rather ad-hoc. Therefore, on physical grounds, we prefer our standard model over the no-scattering case, even though goodness-of-fit 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 CoRoT-1b 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 free-T scenario (see Table D.1) translate, at 95% confidence, to 0.03 <A_B< 0.82, or, when using the best-fit 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 CoRoT-1b (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 model-independent result. Especially, nightside emission measurements are missing, except for our optical phase curve analysis (resulting in a non-detection since the nightside emission is consistent with zero at the level of the measurement errors).

5.3. TrES-2b

Figure 11 shows the various best-fit 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 best-fit parameters as well as 95% credibility regions for the parameters. Parameter posterior distributions are shown in Figs. D.3–D.5.

Results presented in Fig. 12 confirm the very dark nature of TrES-2b, 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).

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 CoRoT-1b. 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).

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.

It is clear that our model provides a relatively good fit to the observed phase curve (Fig. 11, ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}\mathrm{\approx }\mathrm{2}$–2.2, see Table D.2). Compared to the best-fit model by Barclay et al. (2012), our symmetric models consistently calculate a noticeably higher photometric contrast post-eclipse. 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 TrES-2b.

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 day-night 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 post-eclipse 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 post-eclipse maxima are not generally attributed to clouds (e.g., Demory et al. 2013; Esteves et al. 2015). Another possibility for post-eclipse 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.

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 CoRoT-1b, 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 CoRoT-1b (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 CoRoT-1b, we use the calculated A_S values to put constraints on the overall Bond albedo of TrES-2b (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 CoRoT-1b because the spectral coverage is larger and TrES-2b is a cooler planet.

5.4. HAT-P-7b

Figure 16 shows our best-fit models of the different MCMC scenarios. In contrast to TrES-2b, the phase curve of HAT-P-7b 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 best-fit parameters as well as 95% credibility regions for the parameters. Parameter posterior distributions are shown in Figs. D.6–D.8.

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 TrES-2b, 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.

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 ${\mathit{A}}_{\mathrm{G}}\mathrm{=}\frac{\mathrm{2}}{\mathrm{3}}{\mathit{A}}_{\mathrm{S}}$, 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.

As before, ϵ depends on the choice of the thermal model, as illustrated in Fig. 19. Similar to what has been found for TrES-2b, 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 HAT-P-7b, the distribution of the “standard + both” scenario is closer to the “standard + asy” distribution, suggests that for HAT-P-7b, 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⁵ 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 ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}$ 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 under-estimated. However, the ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}$ is dominated by a few points (7 outliers contribute 50% of the total ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}$), 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 under-estimated. When removing the apparently systematic outliers, the ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}$ 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 TrES-2b. These correlations are much stronger and cleaner in this case, since the signal-to-noise ratio of the phase curve is much better for HAT-P-7b.

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 low-albedo 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 CoRoT-1b. However, because of the good data quality of the HAT-P-7b phase curve, the disagreement is stronger than for CoRoT-1b. 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 well-constrained scattering albedo, which is nearly independent of the thermal and asymmetry model that was chosen (see also Fig. 18).

We use the calculated A_S values and measured IR dayside emission spectrum to constrain A_B for HAT-P-7b (Eqs. (36) and (37)). As for TrES-2b, 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 HAT-P-7b is a rather hot planet, compared to TrES-2b. Near-IR 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). Ground-based 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 dust-covered 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 signal-to-noise ratio, spectral resolution or spectral coverage). Even though recent progress has been astonishing, we do not expect exoplanet data to approach Solar-System 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 well-characterized phase curve of Kepler-7b. 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 lon-gitudinal or latitudinal gradients. In reality, this is unlikelyto be true. Secondary-eclipse 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 inter-preted as a non-uniform temperature distribution. IR phasecurves also clearly show temperature gradients (e.g., Knutsonet al. 2007, 2009;Crossfield et al. 2010). In the hypo-thetical no-recirculation limit (ϵ = 0), non-uniformity effects mightproduce a thermal beaming dominated by the sub-stellar pointthat could potentially enhance planetary emitted radiation (e.g.,Selsis et al. 2011; Schwartz &Cowan 2015). However, given the relatively lowsignal-to-noise 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 CoRoT-1b, 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 CoRoT-1b 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 sub-stellar point.

This model has been used to analyze published phase-curve data of CoRoT-1b, TrES-2b and HAT-P-7b. Results are then compared to an analysis of secondary-eclipse data of these planets by Schwartz & Cowan (2015).

We have shown that for CoRoT-1b and HAT-P-7b, inferred albedo and heat recirculation values from optical phase curves are different compared to previously published results. For TrES-2b, both methods yield similar results.

We find that CoRoT-1b 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 CoRoT-1b by Snellen et al. (2009) considered only thermal emission.

In line with previous studies on the optical phase curve of HAT-P-7b (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 re-analyzing further planets with published optical phase curves and reconciling different observations in the optical and the IR for CoRoT-1b.

---

Acknowledgments

This study has received financial support from the French State in the frame of the “Investments for the future” Programme IdEx Bordeaux, reference ANR-10-IDEX-03-02. 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

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.

$\begin{array}{ccc}& & \\ & & \end{array}$Note that α₀ is 0 at opposition and π at primary transit, contrary to α (see Eq. (2)).

For thermal radiation, the phase function is described by the illuminated fraction L of the planetary disk, i.e., the dayside: $\mathrm{\Phi }\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}\cos \mathrm{(}\mathit{z}\mathrm{\left)}\mathrm{\right)}\mathit{.}$(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 10-day orbit around a Sun-like 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 high-frequency structure in the residuals is due to the spatial discretization of the numerical model and has no effect on physical results.

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 2-day orbit around a Sun-like 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 high-frequency 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 signal-to-noise 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 self-consistently, 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: ${\mathit{k}}_{\mathrm{r}}\mathrm{=}\frac{{\mathit{R}}_{\mathrm{p}}}{{\mathit{R}}_{\mathrm{\ast }}}\mathrm{·}$(B.1)Furthermore, for circular orbits, transit duration and transit shape can be related to the orbital inclination i and the projected star-planet separation k_p in units of stellar radii: ${\mathit{k}}_{\mathrm{p}}\mathrm{=}\frac{\mathit{a}}{{\mathit{R}}_{\mathrm{\ast }}}\mathrm{·}$(B.2)Assuming M_P ≪ M_∗, we can write Kepler’s 3rd law as follows: $\frac{{\mathit{P}}_{\mathrm{orb}}^{\mathrm{2}}}{{\mathit{a}}^{\mathrm{3}}}\mathrm{=}\frac{\mathrm{4}{\mathit{\pi }}^{\mathrm{2}}}{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}}\mathrm{·}$(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 HAT-P-7, 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 high-resolution spectroscopy and Van Eylen et al. (2012) who used asteroseismology (see Tables B.1 and B.2 for a compilation).

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 HAT-P-7b 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 best-fit 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 best-fit models result in similar ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}$ values of 10.8 and 10.97, respectively.

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 over-estimated, 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.

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.

$\begin{array}{ccc}& & \\ & & \end{array}$When fitting the phase curve of HAT-P-7b 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 over-estimated. Hence, our results suggest that Eq. (C.2) as a particular choice of phase function is probably not correct for HAT-P-7b. 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 HAT-P-7b most likely has some form of silicate or iron clouds (see, e.g., comparison of Kepler-7b 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).

Appendix C.2: Rayleigh scattering

To investigate the influence of Rayleigh scattering, we incorporated H₂ and He Rayleigh scattering in the model. These two species are thought to form the major constituents of gas-giant atmospheres.

The Rayleigh scattering cross sections of H₂ and He are calculated as ${\mathit{\sigma }}_{\mathrm{ray}\mathit{,i}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{=}{\left(\frac{{\mathit{\lambda }}_{\mathrm{0}\mathit{,i}}}{\mathit{\lambda }}\right)}^{\mathrm{4}}\mathrm{·}\mathit{,}{\mathit{\sigma }}_{\mathrm{0}\mathit{,i}}$(C.3)where σ_ray,i of species i is given in cm² 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 λ₀ 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: ${\mathit{\tau }}_{\mathrm{ray}\mathit{,j}}\mathrm{=}\sum _{\mathit{k}}{\mathit{\sigma }}_{\mathit{k}}{\mathit{C}}_{\mathit{k,j}}\mathit{,}$(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 ${\mathit{C}}_{\mathit{k,j}}\mathrm{=}{\mathit{c}}_{\mathit{k,j}}\mathrm{·}\frac{{\mathit{P}}_{\mathit{k}}\mathrm{-}{\mathit{P}}_{\mathit{k}\mathrm{+}\mathrm{1}}}{{\mathit{\mu }}_{\mathrm{atm}}{\mathit{g}}_{\mathrm{P}}}\mathit{,}$(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₂-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: ${\mathit{\tau }}_{\mathrm{ray}}\mathrm{=}\sum _{\mathit{j}}{\mathit{\tau }}_{\mathrm{ray}\mathit{,j}}\mathit{.}$(C.6)In Fig. C.2, the used Rayleigh scattering cross sections of H₂ 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₂, 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₂ Rayleigh scattering using Penndorf (1957) data is less good.

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, $\begin{array}{ccc}{\mathit{F}}_{\mathrm{R}}& \mathrm{=}& {\mathit{T}}_{\mathrm{ray}}\mathrm{·}{\mathit{F}}_{\mathit{l}}\\ & & \mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{T}}_{\mathrm{ray}}\mathrm{)}\mathrm{·}{\mathrm{\Phi }}_{\mathrm{ray}}\mathrm{·}\frac{\mathit{\omega }}{\cos {\mathit{z}}_{\mathrm{s}}\mathrm{+}\cos {\mathit{z}}_{\mathrm{o}}}\\ & & \mathrm{·}\cos {\mathit{z}}_{\mathrm{s}}\mathrm{·}\frac{\mathit{S}}{\mathit{r}\mathrm{(}\mathit{t}{\mathrm{)}}^{\mathrm{2}}}\mathrm{·}\cos {\mathit{z}}_{\mathrm{o}}\mathrm{·}\mathrm{\Delta \Omega }{\left(\frac{{\mathit{R}}_{\mathrm{p}}}{\mathit{d}}\right)}^{\mathrm{2}}\mathit{,}\end{array}$(C.7)where ω is the single-scattering albedo (set to unity in the all-scattering, zero-absorption approximation used here), T_ray is the transmission along the optical path calculated as ${\mathit{T}}_{\mathrm{ray}}\mathrm{=}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathrm{ray}}\mathrm{(}\frac{\mathrm{1}}{\cos {\mathit{z}}_{\mathrm{s}}}\mathrm{+}\frac{\mathrm{1}}{\cos {\mathit{z}}_{\mathrm{o}}}\mathrm{)}}\mathit{,}$(C.8)with τ_ray the (zenith) optical depth due to Rayleigh scattering (calculated at the mid-point 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 ${\mathrm{\Phi }}_{\mathrm{ray}}\mathrm{=}\frac{\mathrm{3}}{\mathrm{16}\mathit{\pi }}\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{(}\cos {\mathit{\phi }}_{\mathrm{o}}{\mathrm{)}}^{\mathrm{2}}\mathrm{)}\mathit{,}$(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.

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.

Hence, we would expect that Rayleigh scattering does not play a large role, and we neglect it in our phase-curve studies (equivalent to setting P_S = 0). Therefore, phase curves are calculated with the Lambert approximation.

Appendix D: MCMC results

Appendix E: MCMC convergence diagnostics

All Tables

All Figures