---
author:
- 'P. von Paris'
- 'P. Gratier'
- 'P. Bordé'
- 'F. Selsis'
bibliography:
- 'literatur\_idex.bib'
title: 'Inferring heat recirculation and albedo for exoplanetary atmospheres: Comparing optical phase curves and secondary eclipse data'
---

Introduction
============

![image](phase_orbit){width="250pt"} ![image](phase_illu){width="250pt"} ![image](planet_map){width="250pt"}\

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., @sudarsky2000 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 [@sudarsky2000]. The low optical albedos measured for a few hot Jupiters seem to confirm this (e.g., [@rowe2006; @rowe2008]), whereas high albedos inferred for other planets seem to indicate the potential presence of clouds (e.g., [@quintana2013], [@demory2013inhomogen], [@esteves2013]).

Theoretically, the most basic circulation patterns of hot Jupiters are relatively well understood (e.g., [@showman2002], [@madhusudhan2014], [@heng2015]). 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-10mbar). The incident stellar energy is re-radiated before being circulated to the nightside. Consequently, strong temperature contrasts between dayside and nightside develop (e.g., [@showman2002], [@parmentier2013]).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., [@showman2015]). 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., [@showman2011], [@perez2013], [@showman2015]).

Observed thermal phase curves in the (near-)infrared (IR) of a few hot Jupiters seem to confirm these predictions (e.g., [@knutson2007daynight_189733], [@knutson2009daynight_189733], [@stevenson2014]). Assembling many different observations for a large number of hot Jupiters, @cowan2011hot and @schwartz2015 point out a possibly emerging trend, in line with theoretical predictions (e.g., [@perez2013]). Strongly irradiated planets show very inefficient recirculation, whereas for less irradiated planets, the whole range of recirculation is possible.

@cowan2011hot and @schwartz2015 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 @schwartz2015. 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., [@mazeh2010], [@barclay2012], [@esteves2013; @esteves2015]). In such cases, attributing an observed phase curve asymmetry to an offset of the thermal hotspot is physically inconsistent (e.g., [@esteves2015]). Secondly, a few models (e.g., [@snellen2009]) only treat thermal emission and do not include scattered light. Therefore, inferring constraints on the scattering properties of the planet (as done in [@snellen2009]) is equally inconsistent. Thirdly, models that take both thermal and scattering components into account (e.g., [@mislis2012], [@faigler2013], [@placek2014], [@faigler2015]), 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 \[verification\_model\] and eqs. \[cellflux\] and \[thermalflux\]).

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 Sect. \[formodel\] and Sect. \[invmodel\], we describe the physical forward model, the inverse model and the Markov Chain Monte Carlo (MCMC) approach used in this work. Section \[setup\] describes the model setup and planetary scenarios, Sect. \[results\] presents the results, Sect. \[discuss\] a discussion, and we conclude with Sect. \[summary\]. 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.

Forward model {#formodel}
=============

Geometry
--------

The basic geometry set-up of the star-planet system is shown in Fig. \[illu\]. We adopt a coordinate system where the substellar point is fixed throughout the calculations at 0$^{\circ}$ latitude and 180$^{\circ}$ longitude. Therefore, by definition, local latitude $\vartheta$ is 0$^{\circ}$ at the equator (range -90$^{\circ}<\vartheta<90^{\circ}$) and local longitude $\varphi$ is 0$^{\circ}$ at local midnight (range $0^{\circ}<\varphi<360^{\circ}$). 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., [@arras2010] or [@rauscher2014] for a discussion).

The planetary “surface” is divided into cells of 2.5$^{\circ}$x2.5$^{\circ}$ 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 **n**$\cdot$**s**$\geqslant0$ (i.e., dayside) and **n**$\cdot$**o**$\geqslant0$ (i.e., visible to the observer) are satisfied. The nightside is defined with **n**$\cdot$**s**$\leqslant0$, and correspondingly, the part of the nightside visible by the observer with **n**$\cdot$**s**$\leqslant0$ and **n**$\cdot$**o**$\geqslant0$ , simultaneously.

The subobserver latitude $\theta$ is given by the inclination $i$ of the orbital plane with respect to the observer:

$$\label{obslon}
\theta=90^{\circ}-i,$$

i.e., an edge-on orbit has $i$=$90^{\circ}$, and a face-on orbit has $i$=0$^{\circ}$. The subobserver longitude $\phi$ as a function of time $t$ is given by the orbital phase $\alpha$ ($\alpha=0$ and $t=0$ at primary transit, or equivalently, inferior conjunction, see Fig. \[illu\]):

$$\label{obslat}
\phi(t)=360^{\circ} - \alpha(t),$$

with

$$\label{orbphase}
\alpha(t)=T(t)+\omega_P,$$

where $T$ is the true anomaly and $\omega_P$ the argument of periastron. The true anomaly is calculated from the eccentric anomaly $E$:

$$\label{true}
\tan\left(\frac{T}{2}\right)=\sqrt{\frac{1+e}{1-e}}\tan\left(\frac{E}{2}\right),$$

with $e$ eccentricity. $E$ is determined by numerically solving Kepler’s equation:

$$\label{keplereq}
E-e\sin(E)=M,$$

with $M$ mean anomaly. The solution is obtained with a publicly available Fortran procedure[^1], based on a method described in @meeus1991. The mean anomaly $M$ is given by

$$\label{meananomaly}
M=2\pi\cdot \left(x-\lfloor x \rfloor \right),$$

where $x=\frac{t-t_{\rm{peri}}}{P_{\rm{orb}}}$, $P_{\rm{orb}}$ is the planetary orbital period, $t_{\rm{peri}}$ is the time of periastron passage and $\lfloor x \rfloor$ represents the floor function, i.e., the greatest integer less than or equal to $x$. [^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., [@madhusudhan2012] 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 \[rayappend\] 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 **($\mathbf{n} \cdot \mathbf{o} \geqslant 0 $)**, the flux (in Wm$^{-2}$) received by the observer from this cell**, $F_{\rm{r,o,cell}}$,** is given by

$$\label{cellflux}
F_{\rm{r,o,cell}}=\cos z_s \cdot F_{\rm{\ast,p}} \cdot \frac{A_S}{\pi} \cdot  \cos z_o \cdot \frac{\Delta S}{d^2} ,$$

with $z_s$ the local stellar zenith angle, $z_o$ the local observer zenith angle, $F_{\rm{\ast,p}}$ the stellar flux at the planet’s orbit (in W m$^{-2}$), $\Delta S$ the surface element of the cell on the planet (in sr m$^2$), $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. \[cellflux\], the scattering albedo is related to the geometric albedo $A_G$ in a simple manner:

$$\label{ageo} A_G=\frac{2}{3}A_S.$$

The surface element $\Delta S$ of the cell, as seen from the planet’s center, is calculated as

$$\label{surfelement}
\Delta S=\Delta \Omega \cdot R_p^2,$$

with $R_P$ the planetary radius and $\Delta \Omega$ the solid angle (in sr) of the cell (angular extent 2.5$^{\circ}$x2.5$^{\circ}$, see above).

The stellar flux at the planet’s orbit, $F_{\rm{\ast,p}}$, is calculated as follows:

$$\label{stellarflux}
F_{\rm{\ast,p}}=\pi \left(\frac{R_{\ast}}{r}\right)^2 \int_{\lambda_{\rm{low}}}^{\lambda_{\rm{high}}}I_{\rm{\ast,s}}q_I(\lambda)d\lambda,$$

where $I_{\ast,s}$ is a stellar model intensity (in Wm$^{-2}$sr$^{-1}$$\mu$m$^{-1}$) at the stellar surface, $q_I$ is the instrumental filter function, $\lambda_{\rm{low}}$ and $\lambda_{\rm{high}}$ define the wavelength interval of the bandpass, $r$ the star-planet distance and $R_{\ast}$ is the stellar radius[^3]. The numerical integration for eq. \[stellarflux\] is done with a standard trapezoidal integration scheme using a roughly 1nm spectral resolution.

Stellar intensities $I_{\ast,s}$ are obtained from the ATLAS stellar atmosphere grid [@castelli2004][^4]. Instrumental filter functions $T_l$ are taken, for CoRoT-1, from @snellen2009 and, for TrES-2 and HAT-P-7, the Kepler Handbook[^5].

The time-dependent planet-star distance $r$ is calculated with

$$\label{distance}
r(t)=a\frac{1-e^2}{1+e\cos(T(t))},$$

with $a$ the semi-major axis.

The total reflected flux $F_{\rm{r,o}}$ received by the observer is the sum over all contributing cells:

$$\label{refflux}
F_{\rm{r,o}}=\sum_{ (\mathbf{n} \cdot \mathbf{s} \geqslant 0 )\wedge  (\mathbf{n} \cdot \mathbf{o} \geqslant 0 )} F_{\rm{r,o,cell}}.$$

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_{\rm{night}}$ and $T_{\rm{day}}$, respectively.

A widely used approach (for example, [@snellen2009], [@alonso2009], [@mislis2012], [@schwartz2015]) is to relate these temperatures to the Bond albedo $A_B$ and the efficiency of heat re-distribution $\epsilon$. Temperatures are then calculated, based on @cowan2011hot:

$$\label{daytemp}
 T_{\rm{day}}^4=T_{\ast}^4\frac{R_{\ast}^2}{a^2(1-e^2)^{0.5}}\cdot (1-A_B)(\frac{2}{3}-\frac{5}{12}\epsilon)=T_0^4\cdot (\frac{2}{3}-\frac{5}{12}\epsilon),$$

and

$$\label{nighttemp}
 T_{\rm{night}}^4=T_{\ast}^4\frac{R_{\ast}^2}{a^2(1-e^2)^{0.5}}\cdot (1-A_B)\frac{\epsilon}{4}=T_0^4\cdot \frac{\epsilon}{4},$$

where $T_{\ast}$ is the stellar effective temperature and the additional factor $(1-e^2)^{0.5}$ accounts for the mean flux received over an orbit [@williams2002].

The parameter $\epsilon$ is a re-parameterization of the geometrical redistribution factor $f$. @spiegel2010 define $f$ via the apparent dayside flux $F_d$ (i.e., close to secondary eclipse) and the total planetary luminosity $L_P$:

$$\label{f_def}
 L_P=\frac{1}{f}\cdot \pi R_P^2 \sigma T_{\rm{day}}^4=\frac{1}{f}F_d.$$

For perfect heat recirculation, the entire planet is at a uniform temperature, and thus $f=\frac{1}{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 @spiegel2010 and @cowan2011hot, $f=\frac{2}{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. \[f\_def\] to yield

$$\label{f_re}
f=\frac{F_d}{L_P}=\frac{F_d}{F_S}=\frac{T_{\rm{day}}^4}{T_0^4}.$$

For $\epsilon$ 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. \[daytemp\].

Physically, $\epsilon$ 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 $\epsilon$ are expected. For less irradiated planets, a range of circulation regimes is possible (e.g., [@showman2002], [@showman2011], [@perez2013], [@heng2015], [@showman2015]).

Another option in the model is to retain $T_{\rm{night}}$ and $T_{\rm{day}}$ as free parameters for the inverse modeling, an approach used by, e.g., @placek2014.

Since blackbody radiation is isotropic, the thermal flux $F_{\rm{t,o,cell}}$ received by the observer from each cell is given by

$$\label{thermalflux}
F_{\rm{t,o,cell}}(T_c)=\int_{\lambda_{\rm{low}}}^{\lambda_{\rm{high}}}B(T_c,\lambda)q_I(\lambda)d\lambda  \cdot \cos z_o \cdot  \frac{\Delta S}{d^2},$$

with $B(T_c,\lambda)$ the blackbody intensity at the cell’s temperature $T_c$ and $T_c=T_{\rm{night}}$ or $T_c=T_{\rm{day}}$, depending on the location of the cell.

Again, the total emitted flux $F_{\rm{t,o}}$ is the sum over all cells which are visible for the observer (i.e., **n**$\cdot$**o**$\geqslant0$):

$$\label{emmflux}
F_{\rm{t,o}}=\sum_{ (\mathbf{n} \cdot \mathbf{s} \geqslant 0 )\wedge  (\mathbf{n} \cdot \mathbf{o} \geqslant 0 )} F_{\rm{t,o,cell}}(T_{\rm{day}})+\sum_{ (\mathbf{n} \cdot \mathbf{s} \leqslant 0 )\wedge  (\mathbf{n} \cdot \mathbf{o} \geqslant 0 )} F_{\rm{t,o,cell}}(T_{\rm{night}}).$$

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 $\cos z$ in eq. \[cellflux\] compared to eq. \[thermalflux\]. Therefore it is potentially possible to disentangle reflected from emitted light. This effect is illustrated in Fig. \[planckeffect\].

![Effect of thermal emission on the phase curve. See text for discussion.[]{data-label="planckeffect"}](show "fig:"){width="250pt"}\

Figure \[planckeffect\] shows the reflective and the thermal contribution for a hot Jupiter planet in a 2-day orbit around a Sun-like star (500-1,000nm bandpass, $A_B$=$A_S$=0.15, $\epsilon$=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 ($\cos z_s$ in eq. \[cellflux\]), the reflected phase curve is steeper than the thermal one and shows a more pronounced peak towards phase $\alpha=\pi$. 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., [@esteves2013; @esteves2015], [@schwartz2015]) 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, [@agundez2012]) will complicate the interpretation of thermal phase curves. However, because of the currently limited data, assuming uniform hemispheres is justifiable.

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., [@pfahl2008]) and the Doppler boosting (e.g., [@loeb2003]). 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 @quintana2013 to calculate the respective contrasts:

$$\label{ellips}
C_{\rm{ell}}=-A_{\rm{ell}}\cdot \cos(2 \alpha),$$

$$\label{doppler}
C_{\rm{dopp}}=A_{\rm{dopp}}\cdot \sin(\alpha),$$

where $A_{\rm{ell}}$, $A_{\rm{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., [@barclay2012]). 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 [@esteves2013]). $A_{\rm{ell}}$ is given by (see also eq. 2 in [@quintana2013]):

$$\label{ellipsamp}
A_{\rm{ell}}=\alpha_{\rm{ell}}\frac{M_p}{M_{\ast}} \left(\frac{R_{\ast}}{a}\right)^3\sin^2(i),$$

with $M_{\ast}$, $M_p$ the stellar and planetary mass, respectively. $\alpha_{\rm{ell}}$ is a parameter determined by the stellar limb ($u$) and gravity ($g$) darkening:

$$\label{darkening}
\alpha_{\rm{ell}}=0.15\cdot \frac{(15+u)(1+g)}{3-u}.$$

The coefficients $u$ and $g$ depend on stellar characteristics such as effective temperature $T_{\ast}$, metallicity $\rm{[Fe}/\rm{H}]$ and surface gravity $\log g_s$ (determined from $M_{\ast}$ and $R_{\ast}$). As in @barclay2012, @quintana2013 or @esteves2013, we use pre-calculated tables for $u$ and $g$ to interpolate linearly in $T_{\ast}$, $\rm{[Fe}/\rm{H}]$ and $\log g_s$. These tables are taken from model calculations presented in @claret2011 for the Kepler and CoRoT bandpasses. We adopt their coefficients obtained with a microturbulence velocity of 2kms$^{-1}$ with ATLAS stellar models and, in the case of $u$, fitted with a linear least squares approach.

$A_{\rm{dopp}}$ is given by (see also eq. 5 in [@quintana2013]):

$$\label{doppleramp}
A_{\rm{dopp}}=(3-\alpha_{\rm{dopp}})\frac{K}{c},$$

where $c$ is the speed of light, $K$ the radial velocity semi-amplitude and $\alpha_{\rm{dopp}}$ a parameter which depends on the wavelength of observation $\lambda_{\rm{obs}}$ and the stellar effective temperature. As in @quintana2013, we use the approximate equations from @loeb2003 for $\alpha_{\rm{dopp}}$:

$$\label{alphaloeb}
\alpha_{\rm{dopp}}=\frac{e^x(3-x)-3}{e^x-1},$$

where $x=\frac{h\frac{c}{\lambda_{\rm{obs}}}}{k_BT_{\ast}}$ ($h$ Planck’s constant, $k_B$ Boltzmann’s constant). This approach of calculating $(3-\alpha_{\rm{dopp}})$ is somewhat different from the approach used in, e.g., @esteves2013 [@esteves2015]. 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.

$K$ is determined as (with $G$ the gravitational constant)

$$\label{velocityamp}
K=\left(\frac{2\pi G}{P_{\rm{orb}}} \right)^{1/3} \frac{M_p}{M_{\ast}^{2/3}} \sin(i)\frac{1+e\cos \omega_P}{(1-e^2)^{0.5}}.$$

Asymmetries
-----------

As demonstrated by, e.g., @demory2013inhomogen and @esteves2015, 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 re-radiated. Most 3D atmospheric models of hot Jupiters produce such an offset in thermal emission (e.g., [@showman2002], [@showman2015]), and IR phase curves seem to confirm the theoretical predictions (e.g., [@knutson2007daynight_189733]).

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 @demory2013inhomogen. Part of the planet has a scattering albedo $A_S$ (see eq. \[cellflux\]), and part of the planet has a different scattering albedo $d_S\cdot A_S$, with $d_S$ being a free parameter such that $d_S\cdot A_S\leq1$. The extent in longitude of the latter part of the planet is controlled by two further parameters, namely $l_{\rm{start}}$ and $l_{\rm{end}}$. We do not consider any latitudinal variation of the scattering properties.

![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.[]{data-label="asym"}](asym_geometry_scatt "fig:"){width="120pt"} ![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.[]{data-label="asym"}](asym_geometry_therm "fig:"){width="120pt"}\

Second, for the thermal component of the phase curve, we consider a simple offset $\Theta_{\rm{d}}$ of the dayside such that the dayside has an extent in longitude between $\Theta_{\rm{d}}$ and $180^{\circ}+\Theta_{\rm{d}}$. Figure \[asym\] illustrates the asymmetric models.

Phase curve
-----------

The final time-dependent contrast $C$ between star and planet is then

$$\label{contrastref}
C(t)=\frac{F_R(t)+F_E(t)}{F_{\rm{\ast,o}}}+C_{\rm{ell}}+C_{\rm{dopp}}.$$

The stellar flux at the observer, $F_{\rm{\ast,o}}$, is calculated analogous to $F_{\rm{\ast,p}}$ (see eq. \[stellarflux\]):

$$\label{staratobserver}
F_{\rm{\ast,o}}=\pi \left(\frac{R_{\ast}}{d}\right)^2 \int_{\lambda_{\rm{low}}}^{\lambda_{\rm{high}}}I_{\rm{\ast,s}}q_I(\lambda)d\lambda.$$

Model parameter summary
-----------------------

In total, our full physical model as described above contains up to 19 parameters, listed in Table \[parasummary\].

  Group             Parameters
  ----------------- ------------------------------------------------------------
  Stellar (4)       $T_{\ast}$, $R_{\ast}$, $M_{\ast}$, $\rm{[Fe}/\rm{H}]$
  Orbital (4)       $P_{\rm{orb}}$, $e$, $\omega_p$, $i$
  Planetary (2)     $R_p$, $M_p$
  Atmospheric (5)   $A_S$, $A_B$, $\epsilon$, $T_d$, $T_n$
  Asymmetries (4)   $d_S$, $l_{\rm{start}}$, $l_{\rm{end}}$, $\Theta_{\rm{d}}$

  : Parameters of the forward model[]{data-label="parasummary"}

These 19 parameters however are not independent. For instance, if $A_B$ and $\epsilon$ are fixed, $T_d$ and $T_n$ can be determined using eqs. \[daytemp\] and \[nighttemp\].

Inverse model {#invmodel}
=============

  scenario          parameters                                                                            priors                                                                 planets  comments/constraints
  ----------------- ------------------------------------------------------------------------------------- --------------------------------------------------------------------- --------- ------------------------------------------------------
  standard          $\epsilon$, $A_S$, $M_P$                                                              $\epsilon$, $A_S$uniform in \[0,1\]                                    C, T, H  $A_B$=$A_S$
                                                                                                          $M_P$ uniform in \[0.1,10\] $M_{\rm{jup}}$                                      $M_P$ fixed for C
  standard + asy    $\epsilon$, $A_S$, $M_P$,$d_S$, $l_{\rm{start}}$, $l_{\rm{end}}$                      $\epsilon$, $A_S$uniform in \[0,1\]                                     T, H    $A_B$=$A_S$
                                                                                                          $M_P$uniform in \[0.1,10\] $M_{\rm{jup}}$                                       
                                                                                                          $d_S$ uniform in \[0,10\]                                                       $d_S\cdot A_S\leq 1$
                                                                                                          $l_{\rm{start}}$, $l_{\rm{end}}$uniform in \[$90$,$270$\]$^{\circ}$             $l_{\rm{start}} \leq l_{\rm{end}}$
  standard + off    $\epsilon$, $A_S$, $M_P$, $\Theta_{\rm{d}}$                                           $\epsilon$, $A_S$uniform in \[0,1\]                                     T, H    $A_B$=$A_S$
                                                                                                          $M_P$ uniform in \[0.1,10\] $M_{\rm{jup}}$                                      
                                                                                                          $\Theta_{\rm{d}}$uniform in \[0,$360$\]$^{\circ}$                               
  standard + both   $\epsilon$, $A_S$, $M_P$, $\Theta_{\rm{d}}$,$d_S$, $l_{\rm{start}}$, $l_{\rm{end}}$   $\epsilon$, $A_S$uniform in \[0,1\]                                     T, H    $A_B$=$A_S$
                                                                                                          $M_P$uniform in \[0.1,10\] $M_{\rm{jup}}$                                       
                                                                                                          $\Theta_{\rm{d}}$uniform in \[0,$360$\]$^{\circ}$                               
                                                                                                          $d_S$uniform in \[0,10\]                                                        $d_S\cdot A_S\leq 1$
                                                                                                          $l_{\rm{start}}$, $l_{\rm{end}}$uniform in \[$90$,$270$\]$^{\circ}$             $l_{\rm{start}} \leq l_{\rm{end}}$
  free A            $\epsilon$, $A_B$, $A_S$, $M_P$                                                       $A_S$, $A_B$, $\epsilon$ uniform in \[0,1\]                            C, T, H  $a_V \cdot A_S$&lt;$A_B$&lt;$a_V\cdot A_S$+(1-$a_V$)
                                                                                                          $M_P$uniform in \[0.1,10\] $M_{\rm{jup}}$                                       $M_P$ fixed for C
  free T            $T_d$, $T_n$, $A_S$, $M_P$                                                            $A_S$uniform in \[0,1\]                                                C, T, H  
                                                                                                          $T_d$, $T_n$uniform in \[500, 3000\] K                                          $T_n \leq$$T_d$
                                                                                                          $M_P$uniform in \[0.1,10\] $M_{\rm{jup}}$                                       $M_P$ fixed for C
  no scattering     $\epsilon$, $A_B$                                                                     $A_B$,$\epsilon$uniform in \[0,1\]                                        C     $A_S$=0, @snellen2009

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.

$$\label{bayes}
p(V_P | D)\propto p(D | V_P) \cdot p(V_P).$$

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).

MCMC algorithm
--------------

To sample the full parameter space, we adopt a Markov Chain Monte Carlo (MCMC) approach. In this work, we use the emcee python package developed by @foreman2013, implementing an algorithm described in @goodman2010. emcee uses multiple chains (in this work, 500-1,000) 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., [@hou2012]).

To ensure good convergence and avoid any contamination by initial conditions, the chains were run for 500-2,000 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.

![Trace plots of model parameters for the HAT-P-7b “standard + asy” model (see Table \[corot1mcmcsummary\]). Blue line traces the ensemble median, red lines correspond to the \[0.16,0.84\] percentiles (the dark red region in Fig. \[uncertain\_illu\] below).[]{data-label="asym_conv"}](asym_convergence "fig:"){width="250pt"}\

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$\sigma$ uncertainties, respectively. Figure \[uncertain\_illu\] 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 \[corot1mcmcresults\]-\[hat7mcmcresults\], Figs. \[corot1triangle\_2\]-\[hat7triangle\_3\]).

![Example of determination of uncertainty ranges via the cumulative probability distribution: Scattering albedo of CoRoT-1b 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.[]{data-label="uncertain_illu"}](cumulative_alb "fig:"){width="250pt"}\

Goodness-of-fit criteria and model comparison
---------------------------------------------

We use three standard quantities (e.g., [@feigelson2012]) to evaluate the best-fitting models obtained from the MCMC model: $\chi^2$, the reduced $\chi^2_{\rm{red}}$ and the Bayesian Information Criterion (BIC). These are evaluated for the MAP parameter set (not the median of the CDF).

$\chi^2$, the sum of the weighted squared residuals, is defined as

$$\label{chi2def}
\chi^2=\sum_{\rm{j=1}}^{N_D}\left(\frac{M_{V_P,j}-D_j}{\sigma_j}\right)^2,$$

where $N_D$ is the number of data points, $\sigma$ the corresponding uncertainties and $M_{V_P}$ is the model prediction for the parameter vector $V_P$.

The reduced $\chi^2_{\rm{red}}$ is generally calculated by

$$\label{chi2reddef}
\chi^2_{\rm{red}}=\frac{\chi^2}{N_D-N_P},$$

where $N_P$ is the number of parameters. In order to be considered a good fit, $\chi^2_{\rm{red}}$ should be of the order of unity.

Note that by adding more parameters to the fitting model, the $\chi^2$ will decrease. Therefore, we use another criterion, the BIC (valid when $N_D\gg N_P$, which is the case for our calculations), defined as

$$\label{bicdef}
\rm{BIC}=\chi^2 + N_P \ln (N_D).$$

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

$$\label{bicdiff}
p_M=e^{-\frac{\Delta \rm{BIC}}{2}},$$

where $\Delta \rm{BIC}=\rm{BIC}_{M1}-\rm{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_{\rm{model}}$ of a given model:

$$\label{evidencedef}
E_{\rm{model}}=\int p(D | V_P) \cdot p(V_P) dV_P.$$

The calculation of the integral in eq. \[evidencedef\] 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 @goggans2004. However, the calculation is very time-consuming (up to $\sim$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. \[bicdiff\].

Model set-up {#setup}
============

Planets
-------

### CoRoT-1b

CoRoT-1b is the first planet discovered from space [@barge2008] and the first planet with a detected optical phase curve [@snellen2009]. Secondary eclipses of CoRoT-1b have been detected in the optical [@alonso2009] and IR, both from ground (e.g., [@rogers2009] and [@gillon2009]) and from space (e.g., [@deming2011]). These studies suggested that CoRoT-1b would be a very dark planet, with a low albedo ($A_G \lesssim$0.1) and inefficient heat redistribution ($\epsilon \lesssim$0.2), leading to high dayside temperatures of the order of 2,300K.

We use the binned CoRoT red-channel phase curve data presented by @snellen2009 (their Figure 1). All model parameters (stellar parameters, orbital period and inclination, planetary mass and radius) are taken from @barge2008, similarly to what was done in @snellen2009. The orbit is assumed to be circular. This is consistent with the analysis of secondary-eclipse timing (e.g., [@rogers2009], [@alonso2009], [@deming2011]) that constrain the eccentricity to values of $e \lesssim$0.03. Given the data quality of the optical phase curve, the results are not sensitive to eccentricity. When allowing $e$ and $\omega_P$ to be fitted (simulations that are not shown here), constraints on either $A_B$ or $\epsilon$ did not change appreciably compared to the circular case.

### TrES-2b

TrES-2b is a transiting hot Jupiter [@odonovan2006] 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 1,500K (e.g., [@odonovan2010], [@croll2010]). The optical phase curve has been detected in the Kepler data (e.g., [@kipping2011], [@barclay2012], [@esteves2013]). Results indicate a very low geometric albedo ($A_G \leq 0.02$) and a rather efficient day-night redistribution of energy ($\epsilon$&gt;0.5), since dayside and nightside temperatures are quite similar (around 1,300-1,500K, [@esteves2013]).

We use the binned Kepler phase curve data ([@barclay2012], their Figure 5). All model parameters (stellar parameters, orbital period and inclination, planetary radius) are taken from @barclay2012. The orbit is assumed to be circular, since secondary-eclipse timing provided constraints consistent with zero eccentricity (e.g., [@odonovan2010], [@croll2010]) and optical phase curve analysis by @esteves2015 did not find any significant eccentricity.

### HAT-P-7b

As TrES-2b, HAT-P-7b is a very hot Jupiter discovered in the Kepler field [@pal2008] prior to the launch of the satellite. It was the first Kepler planet with a measured optical phase curve [@borucki2009]. Further analysis of the phase curve demonstrated a detection of ellipsoidal variations [@welsh2010]. @christiansen2010 used Spitzer secondary eclipse data to infer maximum brightness temperatures of more than 3,000K. @esteves2013 presented a new optical phase curve using full 3-year Kepler photometry. Measured phase curve amplitudes and secondary eclipse depths vary considerably between @borucki2009, @welsh2010 and @esteves2013, leading to differences in inferred dayside and nightside temperatures of the order of 500K which are most likely a consequence of the extended data set in @esteves2013. Similar brightness temperature differences of 300-500K for the Spitzer secondary eclipse data have been found by @cowan2011hot, compared to @christiansen2010, which they attribute to the use of different stellar models. Based on calculated brightness temperatures and optical phase curves, @christiansen2010 inferred a very inefficient heat redistribution ($\epsilon$ close to zero) between dayside and nightside and modest geometric albedos ($A_G$&lt;0.1). @esteves2013 found a geometric albedo of $A_G=0.18$ and a relatively homogenous temperature distribution, in contrast to @christiansen2010. @schwartz2015 found moderate albedos, slightly lower than @esteves2013 and moderate recirculation values, also in contrast to previous analysis [@christiansen2010].

We use the binned Kepler phase curve data ([@esteves2013], their Figure 3) . However, model parameters are not taken from @esteves2013 or @esteves2015. Instead, stellar parameters are taken from @eylen2012, and planetary parameters (inclination, radius) are taken from @eylen2013. We refer to Appendix \[phase\_transit\] 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., [@pal2008], [@winn2009]). Also, previous phase curve analysis did not find a hint of significant eccentricity [@esteves2015]

Photometric fits
----------------

@snellen2009 analyze the phase curve of CoRoT-1b in terms of the normalized flux $F_{\rm{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:

$$\label{corot1phase}
F_{\rm{norm}}(\alpha)=\frac{1+C(\alpha)}{1+C(0)}.$$

Both @barclay2012 and @esteves2013 fit the phase curves in terms of variations in the photometric light curve, as in eq. \[contrastref\]. They however introduce another, non-physical parameter, a zero-point flux offset $f_0$. 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:

$$\label{f0phase}
L_C(\alpha)=C(\alpha)+f_0.$$

We will use this equation for our analysis of TrES-2b and HAT-P-7b. Hence, $f_0$ is added to the tally of free parameters of the physical model (Table \[corot1mcmcsummary\]).

Planetary scenarios
-------------------

As summarized in Table \[corot1mcmcsummary\], 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. \[daytemp\] and \[nighttemp\] 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., [@snellen2009], [@schwartz2015]) and is consistent with previous studies (e.g., [@alonso2009]). 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\approx0.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:

$$\label{abondconstraint}
a_V \cdot A_S\leq A_B\leq a_V\cdot A_S+(1-a_V).$$

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,\rm{low}}=a_V \cdot 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,\rm{high}}=a_V \cdot 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 \leq T_d$.

To compare directly to the phase-curve analysis by @snellen2009, 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 @snellen2009 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\[corot1mcmcsummary\]), either with asymmetric scattering, a dayside offset or a combination of both.

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. \[formodel\]). The Bond albedo determines the overall received flux which must then be re-emitted on the day and night hemispheres.

$$\label{fluxbalance}
(1-A_B) \cdot F_{\rm{\ast,p}}=2\left (\sum_j^{N_{\rm{day}}}F_{j,\rm{day}}+\sum_j^{N_{\rm{day}}}F_{j,\rm{night}}\right ),$$

where the sums contain the measured fluxes in the various wavelength bands (${N_{\rm{day}}}$ bands of the dayside spectrum, ${N_{\rm{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.

Results
=======

Convergence results
-------------------

The results of the convergence tests and the trace plots for all parameters are shown in Appendix \[convergence\_res\]. 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. \[tres2\_conv\_2\] and \[tres2\_conv\_3\], Table \[tres2gelman\]). Since the parameters describing the albedo asymmetry are tightly anti-correlated, the MCMC algorithm finds basically two separated solutions (see also Fig. \[tres2triangle\_3\]), 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).

CoRoT-1b
--------

Figure \[corot1b\_phase\] shows the data and the best-fit model from @snellen2009 as well as our best-fit models from the various planetary scenarios (Table \[corot1mcmcsummary\]). 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 @snellen2009 provide reasonable fits to the data.

![CoRoT-1b red-channel phase curve: Comparison of best-fit models with data (red) and fit by @snellen2009 (yellow). Primary transit and secondary eclipse not shown.[]{data-label="corot1b_phase"}](corot1_compare "fig:"){width="250pt"}\

Best-fit parameters as well as goodness-of-fit criteria and MCMC posterior parameter distributions are reported in Appendix \[mcmcresults\] (Table \[corot1mcmcresults\], Figs. \[corot1triangle\_2\] and \[corot1triangle\_3\]). Based on the BIC value, the data seem to slightly favor the standard scenario, however, all models in Table \[corot1mcmcresults\] 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 \[corot1mcmcresults\] 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&lt;$A_S$&lt;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. \[corot1b\_scatt\].

![Marginalized posterior distributions for CoRoT-1b scattering albedo in different models. $A_S$ is approximately independent of the thermal model (0.11&lt;$A_S$&lt;0.3 at 95% confidence in the “free A” model).[]{data-label="corot1b_scatt"}](corot1_scattering_albedo "fig:"){width="250pt"}\

Constraints for Bond albedo are weak, and heat recirculation is essentially unconstrained. In addition, constraints on $\epsilon$ show some dependence on the choice of the thermal model employed (see Fig. \[corot1b\_bond\]). 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 [@schlawin2014] which can be interpreted as a hint for the presence of clouds (or, at least, a reflecting layer in the upper atmosphere).

![Marginalized posterior distributions for CoRoT-1b $\epsilon$ (left) and $A_B$ (right) in different models. Constraints on $A_B$ are weak. $\epsilon$ is unconstrained and depends on the thermal model.[]{data-label="corot1b_bond"}](corot1_epsilon "fig:"){width="120pt"} ![Marginalized posterior distributions for CoRoT-1b $\epsilon$ (left) and $A_B$ (right) in different models. Constraints on $A_B$ are weak. $\epsilon$ is unconstrained and depends on the thermal model.[]{data-label="corot1b_bond"}](corot1_bond_albedo "fig:"){width="120pt"}\

Figure \[corot1b\_temp\] 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 @snellen2009, we find that the nightside emission is consistent with zero.

![Marginalized posterior distributions for CoRoT-1b dayside and nightside temperatures in the “free T” model.[]{data-label="corot1b_temp"}](corot1_temp "fig:"){width="250pt"}\

Our results are somewhat contrary to the results from the IR secondary-eclipse measurements ([@schwartz2015]) and conclusions of the analysis of the optical phase curve by @snellen2009. Both studies used the same formalism as our “standard” scenario (i.e., using eqs. \[daytemp\] and \[nighttemp\]), 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 @snellen2009 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. \[corot1b\_circulation\]. For this, we calculated the joint credibility regions in the two-dimensional $A_G$-$\epsilon$ 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 \[corot1b\_circulation\] shows these credibility regions for both the “standard” and the “no scattering” scenario. It is clear that the 1$\sigma$ uncertainty region of @schwartz2015 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 @snellen2009 and @schwartz2015.

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 $\Delta$BIC is small (see Table \[corot1mcmcresults\]).

![Joint credibility regions of recirculation and geometric albedo (see eq. \[ageo\]) for CoRoT-1b in the “standard” (red dots) and “no scattering” (blue dots) scenarios. Orange contour: 1$\sigma$ uncertainty region in @schwartz2015. Our “standard” model strongly disagrees with previous work.[]{data-label="corot1b_circulation"}](Cowan_compare_Corot "fig:"){width="250pt"}\

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. \[abondconstraint\]. Fit results for the free-T scenario (see Table \[corot1mcmcresults\]) translate, at 95% confidence, to 0.03&lt;$A_B$&lt;0.82, or, when using the best-fit values, 0.06&lt;$A_B$&lt;0.8, with $a_V \approx 0.26$ which is consistent with results from the “free A” scenario.

When using the reported IR brightness temperatures of CoRoT-1b ([@rogers2009], [@deming2011]) and our optical brightness temperatures (${N_{\rm{day}}}$=4, ${N_{\rm{night}}}$=1 in eq. \[fluxbalance\]), we obtain $A_B$&lt;0.85, consistent with results from eq. \[abondconstraint\]. 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).

TrES-2b
-------

Figure \[tres2b\_phase\] 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 \[tres2mcmcresults\], we state best-fit parameters as well as 95% credibility regions for the parameters. Parameter posterior distributions are shown in Figs. \[tres2triangle\_1\]-\[tres2triangle\_3\] in the Appendix.

![TrES-2b phase curve: Comparison of best-fit models with data (red) and fit by @barclay2012 (orange). Primary transit and secondary eclipse not shown.[]{data-label="tres2b_phase"}](tres2b_compare "fig:"){width="250pt"}\

Results presented in Fig. \[tres2b\_alb\] confirm the very dark nature of TrES-2b, with scattering albedos $A_S$&lt;0.03 at 95% credibility, as already inferred by previous authors (e.g., [@kipping2011], [@barclay2012], [@esteves2013; @esteves2015]).

![Marginalized posterior distributions for TrES-2b $A_S$ in different models. TrES-2b is a dark planet in all models ($A_S$&lt;0.03 at 95% confidence).[]{data-label="tres2b_alb"}](tres2_scattering_albedo "fig:"){width="250pt"}\

Furthermore, as shown in Fig. \[tres2b\_eps\], the value of $\epsilon$ 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 $\epsilon$ 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 ($\epsilon$ close to zero). As can be seen in Fig. \[tres2b\_eps\], the distributions for $\epsilon$ 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., $\epsilon$ is probably small).

![Marginalized posterior distributions for TrES-2b $\epsilon$ in different models. $\epsilon$ depends strongly on the chosen model.[]{data-label="tres2b_eps"}](tres2_epsilon "fig:"){width="250pt"}\

Figure \[tres2b\_mass\] shows the marginalized posterior distributions for the inferred planetary mass. Also shown are results of previous phase curve modeling by @barclay2012 and @esteves2013 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.

![Marginalized posterior distributions for TrES-2b $M_P$ in different models. Mass from previous studies (including RV measurements) in gray.[]{data-label="tres2b_mass"}](tres2_mass "fig:"){width="250pt"}\

It is clear that our model provides a relatively good fit to the observed phase curve (Fig. \[tres2b\_phase\], $\chi^2_{\rm{red}}\approx$2-2.2, see Table \[tres2mcmcresults\]). Compared to the best-fit model by @barclay2012, our symmetric models consistently calculate a noticeably higher photometric contrast post-eclipse. However, as already noted by @esteves2015, this is simply because the model of @barclay2012 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 [@faigler2015]). This is also the main reason why beaming and ellipsoidal masses do not agree with each other in @barclay2012 or @esteves2013.

When allowing for asymmetric phase curves, the fit becomes slightly better. In terms of the respective BIC values (see Table \[tres2mcmcresults\]), the “standard + off” scenario is slightly favored, although a $\Delta \rm{BIC} \approx 3.5$ is not enough to detect firmly an asymmetry. Our tentative detection is therefore not in contradiction with conclusions of @esteves2015 who state that symmetric models are favored for TrES-2b.

Figure \[tres2triangle\_2\] (right panel) shows some interesting correlations between $\epsilon$, $A_S$ and $\Theta_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 $\epsilon$ are allowed, but that reaches a maximum. Beyond this maximum, scattered light will dominate the phase curve, and again, $\epsilon$ must decrease to produce a significant thermal asymmetry (i.e., large day-night temperature differences).

Figure \[tres2triangle\_3\] (left panel) illustrates a degeneracy in the “standard + asy” scenario, between $A_S$, $d_S$, $l_{\rm{start}}$ and $l_{\rm{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$&lt;1 and $l_{\rm{start}}$ and $l_{\rm{end}}$ delimiting part of the “morning” side, or low $A_S$ and correspondingly $d_S$&gt;1 and $l_{\rm{start}}$ and $l_{\rm{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., [@demory2013inhomogen], [@esteves2015]). 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.

![Joint credibility regions of recirculation and geometric albedo for TrES-2b in the “standard” (black dots) and “standard + off” (blue dots) scenarios. Green contour: 1$\sigma$ uncertainty region in @schwartz2015. Both this and previous work are consistent with each other.[]{data-label="tres2b_circulation"}](Cowan_compare_Tres "fig:"){width="250pt"}\

Figure \[tres2b\_circulation\] shows the constraints on recirculation and geometric albedo as inferred from our “standard” and “standard + off” scenarios, compared to results by @schwartz2015. As for CoRoT-1b, we use joint credibility regions to illustrate our inferred range in the $A_G$-$\epsilon$ plane. Both our results and the results of @schwartz2015 are consistent with each other and certainly agree better than for CoRoT-1b (see Fig. \[corot1b\_circulation\]). However, note that the 1$\sigma$ region of @schwartz2015 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. \[abondconstraint\] and \[fluxbalance\]). Using $a_V \approx 0.4$, we obtain, based on the optical phase curve, $A_B$&lt;0.6. Putting together the dayside brightness temperature measurements from IRAC and Ks bands ([@odonovan2010], [@croll2010]), we then derive $A_B$&lt;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.

HAT-P-7b
--------

Figure \[hat7b\_phase\] 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 \[hat7mcmcresults\], we state best-fit parameters as well as 95% credibility regions for the parameters. Parameter posterior distributions are shown in Figs. \[hat7triangle\_1\]-\[hat7triangle\_3\] in the Appendix.

![HAT-P-7b phase curve: Comparison of best-fit models with data (red) and fit by @esteves2013 (orange). Primary transit and secondary eclipse not shown.[]{data-label="hat7b_phase"}](hat7b_compare "fig:"){width="250pt"}\

Figure \[hat7\_mass\] shows the marginalized posterior distributions for the inferred planetary mass. Also shown are results of previous phase curve modeling by @esteves2013 and @esteves2015 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 @esteves2013 and @esteves2015 (plain and dashed gray lines in Fig. \[hat7\_mass\], 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 [@esteves2013]). Hence, it is clearly demonstrated that a separate fitting of both amplitudes can potentially lead to incorrect mass estimates.

![Marginalized posterior distributions for HAT-P-7b $M_P$ in different models. Mass from previous studies (including RV measurements) in gray.[]{data-label="hat7_mass"}](hat7_mass "fig:"){width="250pt"}\

Figure \[hat7\_alb\] shows the inferred scattering albedo for the different scenarios. The values are broadly consistent with the values from @esteves2015 who find a geometric albedo of $A_G \approx$ 0.2, close to our values (recall $A_G$=$\frac{2}{3} A_S$, eq. \[ageo\]). The fact that $A_S$ is mostly independent of the specific planetary scenario suggests that the estimated value of $A_S$ (0.26&lt;$A_S$&lt;0.34 at 95% confidence) is robust. The arithmetic mean for the combined scenarios is $A_S$=0.28.

![Marginalized posterior distributions for HAT-P-7b $A_S$ in different models. $A_S$ is mostly independent of the adopted model. 0.26&lt;$A_S$&lt;0.34 at 95% confidence.[]{data-label="hat7_alb"}](hat7_scattering_albedo "fig:"){width="250pt"}\

As before, $\epsilon$ depends on the choice of the thermal model, as illustrated in Fig. \[hat7\_eps\]. Similar to what has been found for TrES-2b, the “standard + off” scenario requires a very different distribution for $\epsilon$ 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 \[hat7mcmcresults\]). 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 ($\Delta$BIC$\approx$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 ($\Delta$ BIC&gt;400), again confirming previous phase curve analysis [@esteves2015].

The $\chi^2_{\rm{red}}$ values in Table \[hat7mcmcresults\] 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 $\chi^2_{\rm{red}}$ is dominated by a few points (7 outliers contribute 50% of the total $\chi^2_{\rm{red}}$), 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 $\chi^2_{\rm{red}}$ is reduced to less than 2. This in turn suggests a relatively good fit.

![Marginalized posterior distributions for HAT-P-7b $\epsilon$ in different models. $\epsilon$ is strongly dependent on the adopted planetary scenario.[]{data-label="hat7_eps"}](hat7_epsilon "fig:"){width="250pt"}\

Figure \[hat7triangle\_2\] (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 \[hat7\_circulation\] shows the constraints on recirculation and geometric albedo as inferred from our scenarios, compared to @schwartz2015. As above, we use joint credibility regions to illustrate our inferred range in the $A_G$-$\epsilon$ plane. Note that the “standard+asy” model results and the low-albedo part of the “standard+both” model overlap (see also Fig. \[hat7triangle\_3\]). It is clearly seen that our results and the results of @schwartz2015 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$\sigma$ uncertainty region of @schwartz2015. 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. \[hat7\_alb\]).

![Joint credibility regions of recirculation and geometric albedo for HAT-P-7b in different scenarios. Green contour: 1$\sigma$ uncertainty region in @schwartz2015. Both our and previous work strongly disagree on the inferred $\epsilon$ and $A_G$ values.[]{data-label="hat7_circulation"}](Cowan_compare_Hat "fig:"){width="250pt"}\

We use the calculated $A_S$ values and measured IR dayside emission spectrum to constrain $A_B$ for HAT-P-7b (eqs. \[abondconstraint\] and \[fluxbalance\]). As for TrES-2b, we have $a_V \approx 0.4$. Therefore, based on the optical phase curve, we find 0.11&lt;$A_B$&lt;0.72. The observed Spitzer spectrum translates into $A_B$&lt;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$\mu$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.

Discussion {#discuss}
==========

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., [@arking1968], [@garcia2014], [@petrova2015]), Titan (e.g., [@rages1983]), Jupiter (e.g., [@tomasko1978], [@smith1984]), Saturn (e.g., [@tomasko1984]) or Uranus ([@rages1991], [@pryor1997]).

Large differences are found in the broadband phase curves of, e.g., Jupiter and Saturn (e.g., [@dyudina2005]) or Mars, Mercury and Venus (e.g., [@mallama2009]). 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., [@webber2015], [@hu2015]) 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., [@snellen2009], [@cowan2011hot], [@schwartz2015]):

-   Day and night hemispheres are assumed to be respectively described by a single, uniform temperature, without any longitudinal or latitudinal gradients. In reality, this is unlikely to be true. Secondary-eclipse mapping of the hot Jupiter HD189733b has already demonstrated that the brightness distribution is far from uniform (e.g., [@dewit2012]). This can be interpreted as a non-uniform temperature distribution. IR phase curves also clearly show temperature gradients (e.g., [@knutson2007daynight_189733; @knutson2009daynight_189733], [@crossfield2010]). In the hypothetical no-recirculation limit ($\epsilon=0$), non-uniformity effects might produce a thermal beaming dominated by the sub-stellar point that could potentially enhance planetary emitted radiation (e.g., [@selsis2011], [@schwartz2015]). However, given the relatively low signal-to-noise ratios, the number of effects contributing to the optical phase curve and the large bandpass of Kepler and CoRoT, the optical phase curve is not expected to be very sensitive to the temperature 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., @barclay2012, 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., [@parmentier2013], [@agundez2014], [@showman2015]). 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 @schwartz2015 are marginally compatible, since their 1$\sigma$ uncertainty regions and our 90% credibility regions slightly overlap. Therefore, a re-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.

Conclusions {#summary}
===========

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 re-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 @schwartz2015.

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&lt;$A_S$&lt;0.3 at 95% confidence, which is in slight contrast with previous analyses, which found $A_S$&lt;0.15 ([@snellen2009], [@schwartz2015]). 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 ([@snellen2009], [@schwartz2015]). These contradictions are mainly because previous optical phase curve analysis of CoRoT-1b by @snellen2009 considered only thermal emission.

In line with previous studies on the optical phase curve of HAT-P-7b (e.g., [@esteves2013; @esteves2015]), we find an appreciable albedo ($A_S$$\approx$0.3), slightly higher than inferred from secondary eclipse data. In contrast to previous studies based on secondary eclipse data [@schwartz2015], 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.

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.

Model verification {#verification_model}
==================

We verify that the implementation of model equations leads to the correct limits for reflected and emitted light.

Equations \[lambert\_formula\] and \[z\_formula\] show the phase function for a standard Lambertian sphere. These equations have been used in most studies of optical phase curves so far.

$$\label{lambert_formula}
 \Phi (z)=\frac{1}{\pi}\left(\sin (z)+(\pi-z)\cos(z)\right).$$

$$\label{z_formula}
 \cos (z) = -\sin(i) \cos (\alpha_0).$$

Note that $\alpha_0$ is 0 at opposition and $\pi$ at primary transit, contrary to $\alpha$ (see eq. \[obslat\]).

For thermal radiation, the phase function is described by the illuminated fraction $L$ of the planetary disk, i.e., the dayside:

$$\label{illu_formula}
\Phi (z) = \frac{1}{2}(1- \cos (z)).$$

![Model test: Reflected flux compared to exact Lambertian sphere.[]{data-label="fluxverif"}](veri_GEO "fig:"){width="250pt"}\

In Fig. \[fluxverif\], we show the difference between the exact formulation (eqs. \[lambert\_formula\] and \[z\_formula\]) 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.

![Model test: Emitted flux compared to exact solution.[]{data-label="planckverif"}](veri_THERMAL "fig:"){width="250pt"}\

Figure \[planckverif\] 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, $\epsilon$=0), in order to get an appreciable signal. The difference at peak amplitude is about 0.01ppm for a total amplitude of $\approx$1.6ppm. 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.

Phase curves of transiting planets {#phase_transit}
==================================

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:

$$\label{radius_ratio}
k_r=\frac{R_p}{R_{\ast}}.$$

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:

$$\label{semi_ratio}
k_p=\frac{a}{R_{\ast}}.$$

Assuming $M_P<<M_{\ast}$, we can write Kepler’s 3rd law as follows:

$$\label{kepler3}
\frac{P_{\rm{orb}}^2}{a^3}=\frac{4\pi^2}{GM_{\ast}}.$$

Since the orbital period $P_{\rm{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 \[kepler3\] yields, in this case, an analytic relation between stellar radius and stellar mass.

Such a relation is shown in Fig. \[hat7b\_cons\] (blue line) for HAT-P-7, using a period of $P_{\rm{orb}}$=2.204days [@pal2008] and $k_p$=4.1512 [@esteves2013]. Also shown are stellar parameters taken from @pal2008 who used high-resolution spectroscopy and @eylen2012 who used asteroseismology (see Tables \[hatstellartable\] and \[hattransittable\] for a compilation).

  study        $M_{\ast}$\[$M_{S}$\]    $R_{\ast}$\[$R_{S}$\]  
  ------------ ----------------------- ----------------------- --
  @pal2008     1 .47$\pm$0.08               1.84$\pm$0.23      
  @eylen2012   1.361$\pm$0.021             1.904$\pm$0.01      

  : Stellar parameters for HAT-P-7. In case of asymmetric uncertainties in the original publication, the larger one is stated.[]{data-label="hatstellartable"}

  study                   $k_r$                 $k_p$
  -------------- ----------------------- -------------------
  @pal2008          0.0763$\pm$0.001        4.35$\pm$0.38
  @esteves2013    0.07749$\pm$0.000013    4.1512$\pm$0.0026
  @eylen2013      0.077462$\pm$0.000034   4.1547$\pm$0.0042

  : Planetary parameters for HAT-P-7b. In case of asymmetric uncertainties in the original publication, the larger one is stated.[]{data-label="hattransittable"}

![Consistency between reported radius and mass determinations for the star HAT-P-7, by using the determined orbital period $P_{\rm{orb}}$=2.204days, $a/R_{\ast}$ values from Table \[hattransittable\], blue line is eq. \[kepler3\][]{data-label="hat7b_cons"}](hat7_consistency "fig:"){width="250pt"}\

It is clearly seen that fixing stellar parameters at $R_{\ast}$=1.84$R_{S}$ and $M_{\ast}$=1.47$M_{S}$, as done by @esteves2013, 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 \[corot1mcmcsummary\], i.e., fitting for mass, albedo ($A_B$=$A_S$) and heat recirculation. Consistent models use the stellar parameters from @eylen2012, i.e., $R_{\ast}$=1.90$R_{S}$ and $M_{\ast}$=1.36$M_{S}$, whereas the inconsistent models use $R_{\ast}$=1.84$R_{S}$ and $M_{\ast}$=1.47$M_{S}$, as done in @esteves2013.

Figure \[hat7b\_inconsistent\] shows the residuals $\Delta C=C_{\rm{con}}-C_{\rm{incon}}$ of the best-fit models. As is clearly seen, both scenarios result in virtually identical fits, which only differ by about 0.2ppm (compared to the roughly 75ppm amplitude of the phase curve). Furthermore, both best-fit models result in similar $\chi^2_{\rm{red}}$ values of 10.8 and 10.97, respectively.

![Residuals between consistent (stellar mass and radius from [@eylen2012]) and inconsistent (stellar mass and radius from [@pal2008]) best-fit models of the HAT-P-7b optical phase curve. []{data-label="hat7b_inconsistent"}](hat7b_inconsistent "fig:"){width="250pt"}\

All parameters except planetary mass are not affected by the choice of stellar parameters. Figure \[tricons\] 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 @eylen2012 since these are consistent with parameters deduced from primary transit analysis.

![Marginalized posterior distributions for HAT-P-7b $M_P$, for different adopted stellar parameters in the standard model. []{data-label="tricons"}](mass_inconsistent){width="250pt"}

Optional phase function choices {#rayappend}
===============================

Empirical Solar System phase functions
--------------------------------------

A few previous studies (e.g., [@collier2002], [@kane2010longperiod]) used an empirically derived phase function instead of the Lambert phase function in eq. \[lambert\_formula\]. This phase function was obtained from a fit to optical observations of Venus and Jupiter.

$$\label{delta_m_emp}
   \Delta m (\alpha_0)=0.09\frac{\alpha_0}{100^{\circ}}+2.39\left (\frac{\alpha_0}{100^{\circ}}\right)^2-0.65\left(\frac{\alpha_0}{100^{\circ}}\right)^3.$$

$$\label{emp_phase}
 \Phi (\alpha_0)=10^{-0.4\cdot \Delta m (\alpha_0)}.$$

When fitting the phase curve of HAT-P-7b with the empirical phase function, we obtain a geometric albedo of about $A_G\approx$0.2, consistent with previous estimates using the Lambertian approximation. However, as shown in Fig. \[hat7b\_emp\], the estimated planetary mass is far larger. Even when changing from a uniform prior to a Gaussian prior based on RV measurements (1.78$\pm$0.08, [@esteves2015]), the fitted mass is greatly over-estimated. Hence, our results suggest that eq. \[emp\_phase\] 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 [@webber2015]). Even for Jupiter and Saturn, cloud properties are thought to be responsible for the difference in observed phase functions (e.g., [@dyudina2005]).

Such an impact of the choice of the phase function on mass estimates has also been discussed by @mislis2012.

![Constraints on geometric albedo and planetary mass, as derived from the standard model using the empirical phase function of eq. \[emp\_phase\].[]{data-label="hat7b_emp"}](hat7_empirical_sym "fig:"){width="250pt"}\

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 gas-giant atmospheres.

The Rayleigh scattering cross sections of H$_2$ and He are calculated as

$$\label{rayleigh}
    \sigma_{\rm{ray,i}}(\lambda)= \left (\frac{\lambda_{\rm{0,i}}}{\lambda}\right)^4 \cdot,
    \sigma_{\rm{0,i}}$$

where $\sigma_{\rm{ray,i}}$ of species $i$ is given in cm$^2$ per molecule, $\lambda$ in $\mu$m and $\lambda_{\rm{0,i}}$ is a reference wavelength where $\sigma_{\rm{0,i}}$ has been measured.

This approach is used in many approximative treatments of Rayleigh scattering (see, e.g., [@lecav2008]). For the values of $\lambda_0$ and $\sigma_{\rm{0,i}}$ in eq. \[rayleigh\], measurements from @shardanand1977 were used in this work, as tabulated in Table \[rayleigh\_shard\].

   Molecule   $\lambda_{\rm{0,i}}$ \[$\mu$m\]   $\sigma_{\rm{0,i}}$ \[cm$^2$\]
  ---------- --------------------------------- --------------------------------
    H$_2$                 0.5145                   1.17 $\times$ 10$^{-27}$
      He                  0.5145                   8.6 $\times$ 10$^{-29}$

  : Rayleigh scattering parameters for use in eq. \[rayleigh\][]{data-label="rayleigh_shard"}

The optical depth in an atmospheric layer $j$ due to Rayleigh scattering, $\tau_{\rm{ray,j}}$, is obtained with the following equation:

$$\label{raylayer}
\tau_{\rm{ray,j}}=\sum_k\sigma_k C_{k,j},$$

where $\sigma_k$, $C_{k,j}$ are the Rayleigh cross section and the column density of species $k$ respectively. We calculate the column density as

$$\label{columndens}
C_{k,j}=c_{k,j} \cdot \frac{P_k-P_{k+1}}{\mu_{\rm{atm}} g_P},$$

where $c_{k,j}$ is the volume mixing ratio of species $k$ in layer $j$, $g_P$ planetary gravity, $\mu_{\rm{atm}}$ the mean molecular weight of the atmosphere ($\approx$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 $\tau_{\rm{ray}}$ for use in eq. \[transmission\] is obtained by summing the optical depths of each layer from the surface to the model lid:

$$\label{raytautotal}
\tau_{\rm{ray}}=\sum_j \tau_{\rm{ray,j}}.$$

In Fig. \[h2\_relative\], the used Rayleigh scattering cross sections of H$_2$ are compared to measurements reported in the literature ([@shardanand1977]) as well as different approximations used in various models (Table II of [@penndorf1957], [@lecav2008]).

For H$_2$, the agreement with measurements is very good, again to within the stated error bars of @shardanand1977. Also, the agreement with the approximation of @lecav2008 is very good. The comparison with the parametrization of H$_2$ Rayleigh scattering using @penndorf1957 data is less good.

![Comparison of H$_2$ Rayleigh scattering cross sections. Relative deviations $\frac{\sigma_{\rm{model}}-\sigma_{\rm{data}}}{\sigma_{\rm{model}}}$ in % between model and data sources (as indicated). Vertical lines show measurement uncertainties. Grey horizontal line indicates 0% deviation.[]{data-label="h2_relative"}](rayleigh_cross_relative_H2 "fig:"){width="250pt"}\

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. \[cellflux\] is modified,

$$\begin{aligned}
\label{raycell}
F_R  &=&   T_{\rm{ray}}\cdot F_l+ \\
  \nonumber& & (1-T_{\rm{ray}}) \cdot \Phi_{\rm{ray}} \cdot \frac{\omega}{\cos z_s+\cos z_o} \\
\nonumber  && \cdot \cos z_s \cdot \frac{S}{r(t)^2}\cdot  \cos z_o \cdot \Delta \Omega \left(\frac{R_p}{d}\right)^2 ,\end{aligned}$$

where $\omega$ is the single-scattering albedo (set to unity in the all-scattering, zero-absorption approximation used here), $T_{\rm{ray}}$ is the transmission along the optical path calculated as

$$\label{transmission}
T_{\rm{ray}}=e^{-\tau_{\rm{ray}}(\frac{1}{\cos z_s}+\frac{1}{\cos z_o})},$$

with $\tau_{\rm{ray}}$ the (zenith) optical depth due to Rayleigh scattering (calculated at the mid-point of the spectral interval considered). The value of $\tau_{\rm{ray}}$ will depend critically on the choice of $P_S$ (see above, eqs. \[raylayer\] and \[raytautotal\]). $\Phi_{\rm{ray}}$ is the phase function of Rayleigh scattering

$$\label{rayphase}
\Phi_{\rm{ray}}=\frac{3}{16\pi}\left(1+(\cos\phi_o)^2\right),$$

with $\phi_o$ the angle between observer and incoming stellar light, i.e., $\cos\phi_o=$**s**$\cdot$**o**.

Figure \[rayeffect\] shows the effect of Rayleigh scattering on the phase curve. Together with the standard Lambertian scattering approximation of eq. \[cellflux\], 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$=1bar, the atmosphere starts to become visible, and at $P_S$=10bar, dominates over the “surface” contribution.

![Effect of Rayleigh scattering on the phase curve, for different values of $P_S$. 1ppm error bar to the left. See text for discussion.[]{data-label="rayeffect"}](max_amp_RAY "fig:"){width="250pt"}\

Atmospheric modeling of hot Jupiters predicts the formation of clouds around the 10$^{-2}$bar layer or even at lower pressures (e.g., [@parmentier2013], [@webber2015]). This implies that the reflecting “surface” is at $P_S$&lt;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 @desert2008, Fig. \[tio\] shows the transmission due to TiO and VO absorption assuming $P_S$=10$^{-4}$bar. Figure \[tio\] suggests that not much radiation is expected to penetrate to levels where Rayleigh scattering becomes important.

![Transmission due to TiO and VO, as a function of wavelength. See text for details.[]{data-label="tio"}](tio.png "fig:"){width="250pt"}\

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.

MCMC results {#mcmcresults}
============

![image](corot1_standard_triangle){width="250pt"} ![image](corot1_snellen_triangle){width="250pt"}\

![image](corot1_free_alb_triangle){width="250pt"} ![image](corot1_free_temp_triangle){width="250pt"}\

  scenario        95% credibility regions        best-fit $V_P$      $\chi^2_{\rm{min}}$   $\chi^2_{\rm{red,min}}$  BIC     $\Delta$ BIC   $p_M$
  --------------- ------------------------------ ------------------ --------------------- ------------------------- ------- -------------- -------
  standard        0.06&lt; $A_S$&lt; 0.28        $A_S$= 0.22                19.36                   1.38            24.91   0              1
                  0.18&lt; $\epsilon$&lt; 0.98   $\epsilon$= 0.86                                                                          
  free T          0.13&lt; $A_S$&lt; 0.31        $A_S$= 0.25                17.90                   1.37            26.22   1.31           0.51
                  806&lt; $T_d$&lt;2169K         $T_d$= 874                                                                                
                  523&lt; $T_n$&lt;1634K         $T_n$= 727K                                                                               
  free albedo     0.11&lt; $A_S$&lt; 0.30        $A_S$= 0.25                17.93                   1.37            26.25   1.34           0.51
                  0.11&lt; $A_S$&lt; 0.76        $A_B$= 0.77                                                                               
                  0.04&lt; $\epsilon$&lt; 0.96   $\epsilon$= 0.77                                                                          
  no scattering   0&lt; $A_B$&lt; 0.18           $A_B$= 0.09                24.40                   1.74            29.95   5.04           0.08
                  0&lt; $\epsilon$&lt; 0.31      $\epsilon$= 0.01                                                                          
  @snellen2009    -                              -                          22.24                   1.58            27.79   2.88           0.23

![image](tres2_standard_triangle){width="250pt"} ![image](tres2_free_alb_triangle){width="250pt"}\

![image](tres2b_free_temp_triangle){width="250pt"} ![image](tres2_standard_off_triangle){width="250pt"}\

![image](tres2_standard_asym_triangle){width="250pt"} ![image](tres2_standard_both_triangle){width="250pt"}\

  scenario          95% credibility regions                         best-fit $V_P$                        $\chi^2_{\rm{min}}$   $\chi^2_{\rm{red,min}}$  BIC      $\Delta$ BIC   $p_M$
  ----------------- ----------------------------------------------- ------------------------------------ --------------------- ------------------------- -------- -------------- --------------------
  standard + off    0 &lt; $A_S$&lt; 0.035                          $A_S$= 0.028                                128.48                   2.06            149.42   0              1
                    0.01&lt; $\epsilon$&lt; 0.82                    $\epsilon$= 0.07                                                                                             
                    0.89&lt; $M_P$ &lt; 1.39 $M_{\rm{jup}}$         $M_P$ = 1.09$M_{\rm{jup}}$                                                                                   
                    116&lt; $\theta_{\rm{day}}$&lt; 259$^{\circ}$   $\theta_{\rm{day}}$= 217$^{\circ}$                                                                           
  standard          0 &lt; $A_S$&lt; 0.026                          $A_S$= 0.022                                136.19                   2.19            152.95   3.53           0.17
                    0.11&lt; $\epsilon$&lt; 0.98                    $\epsilon$= 0.94                                                                                             
                    1.08&lt; $M_P$ &lt; 1.48 $M_{\rm{jup}}$         $M_P$ = 1.31$M_{\rm{jup}}$                                                                                   
  standard + asy    0 &lt; $A_S$&lt; 0.06                           $A_S$= 0.046                                126.65                   2.14            155.97   6.55           0.04
                    0.10&lt; $\epsilon$&lt; 0.98                    $\epsilon$= 0.94                                                                                             
                    0.93&lt; $M_P$ &lt; 1.42 $M_{\rm{jup}}$         $M_P$ = 1.06$M_{\rm{jup}}$                                                                                   
                    0.02&lt; $d_S$&lt; 8.87                         $d_S$= 0.14                                                                                                  
                    91&lt; $l_{\rm{start}}$&lt; 230$^{\circ}$       $l_{\rm{start}}$= 98$^{\circ}$                                                                               
                    126&lt; $l_{\rm{end}}$&lt; 267$^{\circ}$        $l_{\rm{end}}$= 192$^{\circ}$                                                                                
  free A            0 &lt; $A_S$&lt; 0.027                          $A_S$= 0.02                                 135.97                   2.22            156.92   7.5            0.02
                    0.03&lt; $A_B$&lt; 0.68                         $A_B$= 0.68                                                                                                  
                    0.04&lt; $\epsilon$&lt; 0.97                    $\epsilon$= 0.86                                                                                             
                    1.09&lt; $M_P$ &lt; 1.49 $M_{\rm{jup}}$         $M_P$ = 1.31$M_{\rm{jup}}$                                                                                   
  free T            0 &lt; $A_S$&lt; 0.028                          $A_S$= 0.020                                135.98                   2.22            156.93   7.51           0.02
                    741&lt; $T_d$&lt;1854K                          $T_d$= 1100K                                                                                                 
                    519&lt; $T_n$&lt;1681K                          $T_n$= 890K                                                                                                  
                    1.09&lt; $M_P$ &lt; 1.48 $M_{\rm{jup}}$         $M_P$ = 1.29$M_{\rm{jup}}$                                                                                   
  standard + both   0 &lt; $A_S$&lt; 0.085                          $A_S$= 0.018                                125.52                   2.16            159.03   9.61           8$\cdot$ 10$^{-3}$
                    0.01&lt; $\epsilon$&lt; 0.95                    $\epsilon$= 0.05                                                                                             
                    0.84&lt; $M_P$ &lt; 1.39 $M_{\rm{jup}}$         $M_P$ = 1.21$M_{\rm{jup}}$                                                                                   
                    0.03&lt; $d_S$&lt; 6.4                          $d_S$= 3.7                                                                                                   
                    92&lt; $l_{\rm{start}}$&lt; 234$^{\circ}$       $l_{\rm{start}}$= 185$^{\circ}$                                                                              
                    127&lt; $l_{\rm{end}}$&lt; 266$^{\circ}$        $l_{\rm{end}}$= 220$^{\circ}$                                                                                
                    81&lt; $\theta_{\rm{day}}$&lt; 318$^{\circ}$    $\theta_{\rm{day}}$= 253$^{\circ}$                                                                           
  @barclay2012      -                                               -                                           125.87                   2.06            146.82   -2.6           3.66

![image](hat7_standard_consistent_triangle){width="250pt"} ![image](hat7_free_alb_consistent_triangle){width="250pt"}\

![image](hat7_free_temp_triangle){width="250pt"} ![image](hat7_standard_off_consistent_triangle){width="250pt"}\

![image](hat7_standard_asym_consistent_triangle){width="250pt"} ![image](hat7_standard_both_triangle){width="250pt"}\

  scenario          95% credibility regions                         best-fit $V_P$                        $\chi^2_{\rm{min}}$   $\chi^2_{\rm{red,min}}$  BIC      $\Delta$ BIC   $p_M$
  ----------------- ----------------------------------------------- ------------------------------------ --------------------- ------------------------- -------- -------------- ----------------------
  standard + asy    0.26&lt; $A_S$&lt; 0.31                         $A_S$= 0.30                                 217.23                   3.68            246.56   0              1
                    0.37&lt; $\epsilon$&lt; 0.99                    $\epsilon$= 0.93                                                                                             
                    1.69&lt; $M_P$ &lt; 1.81 $M_{\rm{jup}}$         $M_P$ = 1.79$M_{\rm{jup}}$                                                                                   
                    0.003 &lt; $d_S$&lt; 0.27                       $d_S$= 0.002                                                                                                 
                    90&lt; $l_{\rm{start}}$&lt; 102$^{\circ}$       $l_{\rm{start}}$= 93$^{\circ}$                                                                               
                    126&lt; $l_{\rm{end}}$&lt; 137$^{\circ}$        $l_{\rm{end}}$= 126$^{\circ}$                                                                                
  standard + both   0.26&lt; $A_S$&lt; 0.35                         $A_S$= 0.34                                 214.74                   3.70            248.26   1.7            0.42
                    0.19&lt; $\epsilon$&lt; 0.93                    $\epsilon$= 0.34                                                                                             
                    1.64&lt; $M_P$ &lt; 1.84 $M_{\rm{jup}}$         $M_P$ = 1.81$M_{\rm{jup}}$                                                                                   
                    0.01 &lt; $d_S$&lt; 0.66                        $d_S$= 0.005                                                                                                 
                    90&lt; $l_{\rm{start}}$&lt; 122$^{\circ}$       $l_{\rm{start}}$= 91$^{\circ}$                                                                               
                    122&lt; $l_{\rm{end}}$&lt; 148$^{\circ}$        $l_{\rm{end}}$= 121$^{\circ}$                                                                                
                    104&lt; $\theta_{\rm{day}}$&lt; 263$^{\circ}$   $\theta_{\rm{day}}$= 259$^{\circ}$                                                                           
  standard + off    0.21&lt; $A_S$&lt; 0.34                         $A_S$= 0.34                                 247.34                   4.05            268.29   21.73          1.9$\cdot$10$^{-5}$
                    0.02&lt; $\epsilon$&lt; 0.44                    $\epsilon$= 0.007                                                                                            
                    1.53&lt; $M_P$ &lt; 1.76 $M_{\rm{jup}}$         $M_P$ = 1.75$M_{\rm{jup}}$                                                                                   
                    118&lt; $\theta_{\rm{day}}$&lt; 233$^{\circ}$   $\theta_{\rm{day}}$= 234$^{\circ}$                                                                           
  standard          0.27&lt; $A_S$&lt; 0.29                         $A_S$= 0.29                                 669.74                   10.80           686.50   439.94         0
                    0.78&lt; $\epsilon$&lt; 0.99                    $\epsilon$= 0.99                                                                                             
                    1.69&lt; $M_P$ &lt; 1.79 $M_{\rm{jup}}$         $M_P$ = 1.75$M_{\rm{jup}}$                                                                                   
  free A            0.27&lt; $A_S$&lt; 0.29                         $A_S$= 0.29                                 669.64                   10.97           690.59   444.03         0
                    0.18&lt; $A_B$&lt; 0.78                         $A_B$= 0.75                                                                                                  
                    0.14&lt; $\epsilon$&lt; 0.99                    $\epsilon$= 0.96                                                                                             
                    1.69&lt; $M_P$ &lt; 1.79 $M_{\rm{jup}}$         $M_P$ = 1.74$M_{\rm{jup}}$                                                                                   
  free T            0.28&lt; $A_S$&lt; 0.29                         $A_S$= 0.29                                 669.68                   10.97           690.62   444.06         0
                    739&lt; $T_d$&lt;2058K                          $T_d$= 1204K                                                                                                 
                    521&lt; $T_n$&lt; 1834K                         $T_n$= 512K                                                                                                  
                    1.70&lt; $M_P$ &lt; 1.79 $M_{\rm{jup}}$         $M_P$ = 1.75$M_{\rm{jup}}$                                                                                   
  @esteves2013      -                                               -                                           276.41                   4.53            297.36   50.8           9.3$\cdot$10$^{-12}$

MCMC convergence diagnostics {#convergence_res}
============================

![image](c1_standard_convergence){width="250pt"} ![image](c1_snellen_convergence){width="250pt"}\

![image](c1_free_alb_convergence){width="250pt"} ![image](c1_free_temp_convergence){width="250pt"}\

![image](t2_standard_convergence){width="250pt"} ![image](t2_free_alb_convergence){width="250pt"}\

![image](t2_asym_convergence){width="250pt"} ![image](t2_free_temp_convergence){width="250pt"}\

![image](t2_off_convergence){width="250pt"} ![image](t2_both_convergence){width="250pt"}\

![image](standard_convergence){width="250pt"} ![image](free_alb_convergence){width="250pt"}\

![image](asym_convergence){width="250pt"} ![image](free_temp_convergence){width="250pt"}\

![image](off_convergence){width="250pt"} ![image](both_convergence){width="250pt"}\

  scenario         $A_S$   $A_B$   $\epsilon$   $T_d$   $T_n$
  --------------- ------- ------- ------------ ------- -------
  standard         1.026     -       1.036        -       -
  free T           1.037     -         -        1.061   1.055
  free albedo      1.029   1.054     1.059        -       -
  no scattering      -     1.024     1.025        -       -

  scenario       $A_S$   $A_B$   $\epsilon$   $T_d$   $T_n$   $M_P$   $\theta_{\rm{day}}$   $l_{\rm{start}}$   $l_{\rm{end}}$   $d_S$
  ------------- ------- ------- ------------ ------- ------- ------- --------------------- ------------------ ---------------- -------
  standard       1.002     -       1.017        -       -     1.002            -                   -                 -            -
  free T         1.031     -         -        1.059   1.09    1.002            -                   -                 -            -
  free albedo    1.002   1.025     1.025        -       -     1.002            -                   -                 -            -
  asymmetric     1.056     -       1.0952       -       -     1.001            -                 1.534             1.251         1.2
  offset         1.002     -       1.022        -       -     1.001          1.022                 -                 -            -
  both           1.041     -       1.164        -       -     1.001          1.148               1.319             1.193        1.218

  scenario       $A_S$   $A_B$   $\epsilon$   $T_d$   $T_n$   $M_P$   $\theta_{\rm{day}}$   $l_{\rm{start}}$   $l_{\rm{end}}$   $d_S$
  ------------- ------- ------- ------------ ------- ------- ------- --------------------- ------------------ ---------------- -------
  standard       1.006     -       1.013        -       -     1.002            -                   -                 -            -
  free T         1.093     -         -        1.137   1.124   1.002            -                   -                 -            -
  free albedo    1.01    1.061     1.037        -       -     1.001            -                   -                 -            -
  asymmetric     1.025     -       1.112        -       -     1.002            -                 1.009             1.011        1.021
  offset         1.099     -       1.063        -       -     1.002          1.167                 -                 -            -
  both           1.089     -       1.256        -       -     1.002          1.332               1.167             1.257        1.186

[^1]: www.davidgsimpson.com/software/keplersoln\_f90.txt

[^2]: In practice, we first calculate $M$, $T$ and $\alpha$ with $t_{\rm{peri}}=0$ and then interpolate such that $t=0$ occurs at $\alpha=0$.

[^3]: Note that we calculated the star’s solid angle, as seen from the planet, in eq. \[stellarflux\] as $\pi \left(\frac{R_{\ast}}{r}\right)^2$. This assumes that $R_{\ast}\ll r$, a condition that is on the verge of breaking down for close-in 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.

[^4]: http://www.user.oats.inaf.it/castelli/grids.html

[^5]: retrieved from http://keplergo.arc.nasa.gov/Instrumentation.shtml

[^6]: Note that for tidally-locked 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.
