Closed-formed solutions of geometric albedos and phase curves of exoplanets for any reﬂection law

The albedo of a celestial body is the fraction of light reﬂected by it. Studying the albedos of the planets and moons of the Solar System dates back at least a century [1, 2, 3, 4, 5]. Of particular interest is the relationship between the albedo measured at superior conjunction (full phase), known as the “geometric albedo”, and the albedo considered over all phase angles, known as the “spherical albedo” [2, 6, 7]. Modern astronomical facilities enable the measurement of geometric albedos from visible/optical secondary eclipses [8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20] and the inference of the Bond albedo (spherical albedo measured over all wavelengths) from infrared phase curves [21, 22, 23, 24, 25] of transiting exoplanets. Determining the relationship between the geometric and spherical or Bond albedos usually involves complex numerical calculations [26, 27, 28, 29, 30, 31, 32] and closed-form solutions are restricted to simple reﬂection laws [33, 34]. Here we report the discovery of closed-form solutions for the geometric albedo and integral phase function that apply to any law of reﬂection. The integral phase function is used to obtain the phase integral, which is the ratio of the spherical to the geometric albedos. The generality of the solutions stems from a judicious choice of the coordinate system in which to perform diﬀerent parts of the derivation. The closed-formed solutions have profound implications for interpreting observations. The shape of a reﬂected light phase curve and the secondary eclipse depth may now be self-consistently inverted to retrieve fundamental physical parameters (single-scattering albedo, scattering asymmetry factor, cloud cover). Fully-Bayesian phase curve inversions for reﬂectance maps and si-multaneous light curve detrending may now be performed, without the need for binning in time, due to the eﬃciency of computation. We demonstrate these innovations for the hot Jupiter Kepler-7b, in-ferring a revised geometric albedo of 0 . 12 ± 0 . 02, a Bond albedo of 0 . 18 ± 0 . 03 and a phase integral of 1 . 5 ± 0 . 1, which is consistent with isotropic scattering. The scattering asymmetry factor is 0 . 04 ± 0 . 15, implying that the aerosols are small compared to the wavelengths probed by the Kepler space telescope. In the near future, one may use the closed-form solutions discovered here to extract fundamental parameters, across wavelength, from multi-wavelength phase curves of both gas-giant and terrestrial exoplanets measured by the James Webb Space Telescope.

The albedo of a celestial body is the fraction of light reflected by it. Studying the albedos of the planets and moons of the Solar System dates back at least a century [1,2,3,4,5]. Of particular interest is the relationship between the albedo measured at superior conjunction (full phase), known as the "geometric albedo", and the albedo considered over all phase angles, known as the "spherical albedo" [2,6,7]. Modern astronomical facilities enable the measurement of geometric albedos from visible/optical secondary eclipses [8,9,10,11,12,13,14,15,16,17,18,19,20] and the inference of the Bond albedo (spherical albedo measured over all wavelengths) from infrared phase curves [21,22,23,24,25] of transiting exoplanets. Determining the relationship between the geometric and spherical or Bond albedos usually involves complex numerical calculations [26,27,28,29,30,31,32] and closed-form solutions are restricted to simple reflection laws [33,34]. Here we report the discovery of closed-form solutions for the geometric albedo and integral phase function that apply to any law of reflection. The integral phase function is used to obtain the phase integral, which is the ratio of the spherical to the geometric albedos. The generality of the solutions stems from a judicious choice of the coordinate system in which to perform different parts of the derivation. The closed-formed solutions have profound implications for interpreting observations. The shape of a reflected light phase curve and the secondary eclipse depth may now be selfconsistently inverted to retrieve fundamental physical parameters (single-scattering albedo, scattering asymmetry factor, cloud cover). Fully-Bayesian phase curve inversions for reflectance maps and simultaneous light curve detrending may now be performed, without the need for binning in time, due to the efficiency of computation. We demonstrate these innovations for the hot Jupiter Kepler-7b, inferring a revised geometric albedo of 0.12 ± 0.02, a Bond albedo of 0.18 ± 0.03 and a phase integral of 1.5 ± 0.1, which is consistent with isotropic scattering. The scattering asymmetry factor is 0.04±0.15, implying that the aerosols are small compared to the wavelengths probed by the Kepler space telescope. In the near future, one may use the closedform solutions discovered here to extract fundamental parameters, across wavelength, from multiwavelength phase curves of both gas-giant and terrestrial exoplanets measured by the James Webb Space Telescope.
Keywords: exoplanetary atmospheres, albedos, reflected light phase curves, phase integral, analytical formulae, data analysis and interpretation

Main Text
Named after the astronomer George Bond [1], the Bond albedo A B is the fraction of starlight that is reflected by a celestial body-either by its atmosphere or surface-over all viewing angles and wavelengths of light. For example, the Bond albedo of the Earth is about 0.3. The spherical albedo A s phase angle exoplanet α Observer's perspective An observer on Earth views an exoplanet with observer-centric latitude Θ and longitude Φ. At each point in its orbit, the exoplanet resides at a phase angle α. Starlight impinges upon the exoplanet at an angle θ from its zenith. An infinitesimal area element on the exoplanet may be described by the exoplanetcentric coordinate system (θ, φ). Trios of these angles form triangles on the surface of the sphere, and may be related by the spherical law of cosines.
is considered over all angles, but only at one wavelength of light. The geometric albedo A g is the albedo measured when an exoplanet is at superior conjunction or full phase-when the star is between the observer on Earth and the exoplanet. This occurs at an orbital phase angle α of zero. The phase integral q relates the spherical and geometric albedos [2,6], The phase integral is notoriously difficult to evaluate [7] and requires knowledge of the integral phase function Ψ [2,6], An exact, closed-form solution for Ψ exists for a Lambertian sphere [2, 6]-a perfectly reflective sphere that is equally bright regardless of viewing direction. Insight is gained by considering the geometry of the problem [3,6,26], which may be formulated in two different coordinate systems ( Figure 1). The observer-centric coordinate system defines a latitude Θ and longitude Φ, such that the location on the exoplanet closest to Earth (i.e., the sub-Earth point) always sits at Φ = 0. The position of the exoplanet in its orbit around the star is defined by a phase angle α, such that the exoplanet is illuminated only across α − π/2 ≤ Φ ≤ π/2 as seen by the observer. The exoplanet-centric coordinate system defines a polar angle θ and an azimuthal angle φ for each point on the exoplanet. Starlight impinges upon each point at a zenith angle θ . The observer-and exoplanet-centric coordinate systems are generally not aligned in the same plane (Figure 1). In principle, calculations may be performed in either coordinate system and all six angles are mathematically related (see Methods).
All of the quantities of interest (A g , A s , Ψ, q) fundamentally involve the reflection coefficient [6,7,27,28,33], which relates the intensities of incident starlight I and reflected starlight I 0 . We have defined µ = cos θ . By definition, a Lambertian sphere has ρ = 1 [6]. The geometric albedo may be evaluated in the exoplanet-centric coordinate system [6,27], where we have defined µ = cos θ. The reflection coefficient evaluated at zero phase angle (α = 0) is represented by ρ 0 . The integral phase function may be written as Ψ = F/F 0 and is evaluated in the observer-centric coordinate system using [6] F = π/2 α−π/2 π/2 0 ρµµ cos Θ dΘ dΦ, where F 0 is F evaluated at α = 0. Alternatively, the geometric albedo may be evaluated in the exoplanet-centric coordinate system as A g = 2F 0 /π [6]. The reflection coefficient is derived from the radiative transfer equation (see Methods). If the approximation is made that single scattering may be described by any reflection law, but multiple scattering of light occurs isotropically [35], then the reflection coefficient of a semi-infinite atmosphere is where H(µ) is Chandrasekhar's H-function [36] and H = H(µ ) (see Methods). The single-scattering albedo ω is the fraction of light reflected during a single scattering event. The reflection law or scattering phase function P relates the incident direction of collimated starlight to the emergent direc- For illustration, the scattering asymmetry factor of the Henyey-Greenstein reflection law is chosen to be g = 0.508 such that the phase integral (single scattering only) has a value of 3.5, which matches the ensemble averaged value of a population of hot Jupiters [24]. The Lambertian sphere has geometric and spherical albedos of 2/3 and 1, respectively.
tion of scattered starlight (see Methods). It is usually assumed to depend only on the difference in angle between these directions, which is known as the scattering angle β [37]. For example, isotropic and Rayleigh scattering are described by P = 1 and P = 3 4 (1 + cos 2 β), respectively. The scattering and phase angles are trivially related, cos β = − cos α.
At zero phase angle, one simply has cos β = −1. Correspondingly, the scattering phase function at α = 0 is evaluated at µ = µ and φ = π (see Methods) and is represented by P 0 . For example, P 0 = 3/2 for Rayleigh scattering. Since P 0 is a number and not a function, the geometric albedo for any scattering law is straightforwardly derived (see Methods), where = (1 − γ)/(1 + γ) is the bihemispherical reflectance [35] and γ = √ 1 − ω. The bihemispherical reflectance is typically small unless one approaches the pure reflection limit (ω = 1). The term involving P 0 accounts for the single scattering of light, while the other terms account for isotropic multiple scattering. For ω = 1 and isotropic single scattering (P 0 = 1), one obtains A g = 17/24, which is slightly higher than the A g = 2/3 value of the Lambertian sphere by 1/24.
The integral phase function Ψ may be derived using the stated reflection coefficient. In the observercentric coordinate system, the two-dimensional integral over Θ and Φ may be evaluated as two independent, one-dimensional integrals. In the exoplanet-centric coordinate system, α has a complicated relationship with θ, θ and φ. In the observer-centric coordinate system, α is instead the third independent angle, which allows P (α) to be taken out of the integrals involving Θ and Φ. These insights render the derivation of Ψ analytically tractable (see Methods).
Physically, Ψ consists of three components. The first component Ψ S exhibits single-scattering behaviour. The second component Ψ L behaves like a Lambertian sphere [35]. The third component Ψ C is a correction term that is significant only when ω ∼ ∼ 1 (strong reflection). Our derivation yields where the various coefficients are and (b) ω = 1. As the strength of scattering increases, the integral phase function becomes approximately Lambertian. For illustration, the scattering asymmetry factor of the Henyey-Greenstein reflection law is chosen to be g = 0.508 such that the phase integral (single scattering only) has a value of 3.5, which matches the ensemble averaged value of a population of hot Jupiters [24].
The three components of the integral phase function are and we have defined α ± = sin(α/2)±cos(α/2). For −π ≤ α ≤ 0, one needs to replace α by |α|. While these formulae may appear long and unwieldly, they constitute a closed-form solution of the integral phase function for any reflection law P that includes both single and (isotropic) multiple scattering. The phase integral q may be straightforwardly obtained by a numerical integration of Ψ over the phase angle. The spherical albedo is A s = qA g (see Methods for an alternative way of computing A s ). Figure 2 shows calculations of A g , A s and q for the isotropic, Rayleigh and Henyey-Greenstein reflection laws. Multiple scattering is dominant in determining the values of the geometric and spherical albedos. For example, if only isotropic single scattering is considered, the maximum values of the geometric and spherical albedos are 1/8 and 2(1 − ln 2)/3 ≈ 0.2, respectively. If only single scattering is considered, the phase inte-gral is independent of the single-scattering albedo ω. With multiple scattering included, q has a nontrivial dependence on ω, especially for the Henyey-Greenstein reflection law. We validated our calculations of A g and A s against those from previous studies [27,33,34] for various reflection laws (see Methods and Figure 5).
There are profound implications of the closedform solutions of A g and Ψ for fitting and interpreting the measured reflected light phase curves of exoplanets. Let the fluxes of the exoplanet and star observed at Earth be F p and F , respectively. The reflected light phase curve is [7] where R is the radius of the exoplanet and a is its orbital semi-major axis. Figure 3 shows several examples of Ψ. At low and intermediate values of the single-scattering albedo ω, the shapes of reflected light phase curves encode information on the reflection laws. As ω approaches unity, the phase curve shape follows that of a Lambertian sphere [34]. A severe test of the theory is to measure reflected light phase curves of cloudfree exoplanetary atmospheres, where Rayleigh scattering of light is mediated by atoms and molecules and the phase curve is predicted to have a distinct shape if ω < 1 ( Figure 3). We use the closed-form solutions to fit the phase curve of the hot Jupiter Kepler-7b measured by the   : the baseline single-scattering albedo ω0, the total single-scattering albedo ω, the scattering asymmetry factor g, the local longitudes x1 and x2 between which the single-scattering albedo is ω0, the power in the spherical harmonic mode describing the temperature map for thermal emission C11; and posterior distributions for the derived parameters (blue labels), which can be computed from the fundamental parameters: the geometric albedo Ag, the Bond albedo AB, the phase integral q, and the ratio of flux in reflected light at secondary eclipse F0 to the maximum reflected light Fmax.
Kepler space telescope. Since the reflected light component of the phase curve has a peak offset [38], a model of an inhomogeneous atmosphere [39,40] is needed to fit the data (see Methods). For tidallylocked hot Jupiters, an inhomogeneous atmosphere is a natural consequence of aerosols having finite condensation temperatures [41,42]. It is assumed that the atmosphere has a single-scattering albedo ω 0 between two local longitudes x 1 and x 2 , while its other regions are enhanced by ω such that the total albedo is ω = ω 0 + ω ( Figure 6). The asymmetry of scattering is globally quantified by the scattering asymmetry factor g (see Methods). Standard phase curve analyses perform detrending of the raw telescope measurements and fitting to a phase curve model as separate steps [15,16,20,38]. Often after detrending, the photometry is binned in time to reduce data volume and increase the signal-to-noise of the weak phase curve signal. The shape of the phase curve, geometric albedo, and secondary eclipse depth are usually treated as independent fitting parameters. A key improvement from previous analyses is the ability to jointly and self-consistently fit the shape of the phase curve and secondary eclipse depth in terms of fundamental physical parameters (ω 0 , ω, g, x 1 and x 2 ). A major implication is the ability to simultaneously detrend the phase curve in the time domain and fit for physical parameters at the native temporal resolution of the photometry-without the need to bin in time (see Methods). Performing both steps simultaneously and without binning [43] ensures that uncertainties in the detrending process are accurately propagated into the physical parameters. Figure 4 shows the model fit to the Kepler-7b phase curve, which indicate that the atmosphere has a dark region with ω 0 = 0.04 +0.05 −0.03 between the local longitudes of x 1 = −12 +8 −9 and x 2 = 80 +7

−11
degrees. The reflective regions of the atmosphere have a single-scattering albedo of ω = 0.20 +0.07 −0.09 and a scattering asymmetry factor of g = 0.04 +0. 15 −0.13 , implying roughly isotropic scattering of starlight by aerosols with sizes smaller than the wavelengths probed by the Kepler bandpass (0.42-0.90 µm). It is important to note that the inference of g ≈ 0 is not a consequence of assuming isotropic multi-ple scattering as ω is considerably less than unity. The geometric albedo is A g = 0.12 ± 0.02, which significantly revises the previously reported value of 0.35 ± 0.02 [13,38]. Using the inferred values of ω 0 , ω, g, x 1 and x 2 , the Bond albedo and phase integral are predicted to be A B = 0.18 ± 0.03 and q = 1.46 +0.10 −0.08 , respectively. The transition between the dark and reflective regions of the atmosphere occurs at ∼ 2000 K, which we interpret as the condensation temperature of the aerosols. A broad class of aerosol compositions, including the olivine group, are ruled out for solar metallicity [44].
In the upcoming era of the James Webb Space Telescope, phase curves of exoplanets from 0.6 to 24 µm will be procured. The closed-form solutions presented here enable fundamental physical parameters to be extracted from multi-wavelength phase curves. Alternatively, one may construct phase curve models where the single-scattering albedo and scattering asymmetry factor are inserted, from first principles (using Mie theory), as functions of wavelength. Both approaches will allow phase curve models to be rigorously and efficiently confronted by data and open up new avenues of cloud/haze research in the atmospheres of exoplanets.

Radiative transfer equation
Let the scattering phase function be represented by P , which describes the mathematical relationship between the incident and emergent angles of radiation. It is chosen to have no physical units; an alternative choice is for it to have units of per unit solid angle. The fraction of light reflected during a single scattering event is given by the single-scattering albedo ω. With these definitions, the radiative transfer equation reads [6,36,37,45,46,47,48] The optical depth τ is the generalisation of length in radiative transfer and takes into account not just the spatial extent of the atmosphere, but how dense it is and the ability of its constituent atoms, molecules, ions and aerosols to absorb and scatter radiation. Since we are interested in calculating the intensity of reflected light, equation (14) ignores the contribution of thermal (blackbody) emission. The term involving the integral accounts for the multiple scattering of radiation, while the last term describes the collimated beam of incident starlight.
The integral is generally difficult to evaluate, because the integration is performed over all incident angles of radiation (with dΩ denoting the incident solid angle). The scattering phase function is normalised such that [37,47] 1 4π with dΩ denoting the emergent solid angle. It depends on both the incident and emergent angles: P (µ , µ, φ). For the stellar beam, we have P = P (−µ , µ, φ); the minus sign accounts for the convention chosen that radiation travelling into the atmospheric column, towards the center of the exoplanet, corresponds to µ < 0. At zero phase angle, the scattering phase function becomes P 0 = P (−µ , µ , π). The scattering asymmetry factor is the first moment of the scattering phase function [37,47], 2π 0 π 0 P cos β sin β dβ dφ.

Solutions of radiative transfer equation
If one ignores the integral in equation (14), then only single scattering is considered. For a semiinfinite atmosphere, the intensity of reflected light is [35] In the limit of isotropic single scattering (P = 1), one obtains the Lommel-Seeliger law of reflection [49,50]. For isotropic multiple scattering, Chandrasekhar's exact solution ignoring the stellar beam is [3,36] where H = H(µ ) and H(µ) is the Chandrasekhar H-function [36], An identical solution for I 0,M may be found by evaluating the integral in equation (14) using the twostream solutions, but with the H-function taking on an approximate form (Hapke's solution) [35,51], where we have γ = √ 1 − ω. To combine the solutions, one needs to account for isotropic multiple scattering of the stellar beam, which yields [35]

Region of enhanced reflectivity
Region of normal reflectivity α + x 2 α + x 1 Figure 6: Regions of normal and enhanced reflectivity for an inhomogeneous atmosphere in terms of the longitude Φ within the observer-centric coordinate system. In the local longitude of the exoplanet (where the substellar point sits at x = 0), the atmosphere has a single-scattering albedo of ω0 across x1 ≤ x ≤ x2. When this region is within view of the observer, it is highlighted with a thick blue line. Regions of enhanced reflectivity (with a single-scattering albedo of ω = ω0 + ω ) that are within the observer's view are highlighted with thick red lines.
The reflection coefficient is obtained simply using ρ = I 0 /I µ .

Formulae of albedos
If only single scattering is considered, it is straightforward to derive the geometric albedo, The derivation may be performed in either the exoplanet-or observer-centric coordinate systems. To evaluate the geometric albedo for isotropic multiple scattering, one uses Hapke's linear approximation of the Chandrasekhar H-function [35], The geometric albedo associated with isotropic multiple scattering only is The total geometric albedo is A g = A g,S + A g,M .
There is an alternative way of deriving the spherical albedo. In Figure 1, the flux reflected by the infinitesimal area element with projected solid angle µ dΩ is I 0 µ dΩ = ρI µ µ dΩ. In the exoplanetcentric coordinate system, one has dΩ = dµ dφ.
The incident flux of starlight is πµ I . Therefore, the plane albedo associated with each atmospheric column is [6] where 0 ≤ µ ≤ 1 represents only the emergent (and not downwelling) radiation from the exoplanet. The fraction of incident starlight reflected by the entire exoplanet needs to consider all of the atmospheric columns. The reflected flux is 2π µ πI A dµ , while the incident flux is π 2 I . The ratio of these two quantities yields the spherical albedo [6], The isotropic multiple scattering component of the spherical albedo is [35] A s,M = 5 6 + where Hapke's linear approximation of the Chandrasekhar H-function has again been invoked. Finally, the spherical albedo including both single and isotropic multiple scattering is where the single scattering phase integral q S is obtained from For isotropic single scattering (P = P 0 = 1), we verified numerically that q S = 16(1 − ln 2)/3 ≈ 1.63655 [2,28]. For Rayleigh scattering, the singlescattering component of the spherical albedo may be derived exactly, which has a value of about 0.208144ω. We verified numerically that q S ωP 0 /8 ≈ 0.208144ω.
In the limit of conservative (ω = 1) Rayleigh single scattering and isotropic multiple scattering, we obtain which is not inconsistent with the estimate of 0.751 for conservative Rayleigh scattering [27] that is the basis for the commonly quoted factor of 3/4 [29,34]. The slight discrepancy of about 2.3% stems from the assumption of isotropic multiple scattering.

Dimensionless Sobolev fluxes
Using Hapke's linear approximation of the Chandrasekhar H-function [35], the reflection coefficient may be written as where the various coefficients have already been defined in equation (10). The dimensionless flux of Sobolev may be written as three terms, For a homogeneous atmosphere, the term that exhibits single scattering behaviour is The Lambertian-sphere-like term is The correction term is Two key properties of all three terms are worth belabouring: the integrals over Θ and Φ may be evaluated independently, and closed-form solutions exist for them. By defining F 0 = F (α = 0), the integral phase function of a homogeneous atmosphere may be obtained using Ψ = F/F 0 .

Henyey-Greenstein reflection law
The widely used Henyey-Greenstein scattering phase function is [52] where the scattering angle β is [37,47,48] cos The scattering asymmetry factor g quantifies the asymmetry of scattering: g = −1, 0 and 1 correspond to purely reverse, isotropic and purely forward scattering, respectively. The scattering phase function at zero phase angle (cos β = − cos α = −1) is

Validation
To validate our calculations of A g and A s , we compare them against the numerical calculations of Dlugach & Yanovitskij [27], who tabulated the geometric and spherical albedos from ω = 0.7 to 1 ( Figure 5). For isotropic and Rayleigh scattering, the discrepancies associated with A g and A s are typically ∼ 0.1-1%. For the Henyey-Greenstein law of reflection, the discrepancies associated with A g and A s are ∼ 10% for g = 0.5. The larger discrepancies indicate that the assumption of isotropic multiple scattering starts to break down when scattering becomes strongly asymmetric. For isotropic scattering and non-conservative (ω = 1) Rayleigh scattering, we perform further comparisons of the spherical albedo to the empirical fitting functions of van de Hulst [33] and Madhusudhan & Burrows [34], respectively. The differences between these fitting functions and our calculations are negligible (< 1%), as shown in Figure 5.

Inhomogeneous atmosphere
Let the local longitude of the exoplanet be x. Since it lies in the same plane as α, it is trivially related to the phase angle by x = Φ − α [39]. It is assumed that the region from x 1 ≤ x ≤ x 2 has a single-scattering albedo of ω 0 , while all of the other regions have an enhanced single-scattering albedo of ω 0 + ω [39]. The dimensionless Sobolev fluxes are generalised to The coefficients ρ S , ρ L and ρ C are given by equation (10), while Ψ S , Ψ L and Ψ C are given by equation (11). The total dimensionless Sobolev flux is F = F S + F L + F C , while F 0 = F (α = 0) and A g = 2F 0 /π. The integral phase function, Ψ = F/F 0 , is generally not symmetric about α = 0 and the generalised expression for the phase integral does not assume this symmetry [39], To compute Ψ S , Ψ L and Ψ C , one needs to determine the regions of the inhomogeneous atmosphere that are within view of the observer [39] in the observer-centric coordinate system ( Figure 6). If we write Ψ i for short (with the index i = S, L or C), then we can express Ψ i in terms of the following functions,  Figure 6, one has to construct Ψ i according to the correct sextant, which is expressed as a set of inequalities of Φ. In comparison to a homogeneous atmosphere, an inhomogeneous atmosphere has ω , x 1 and x 2 as additional parameters.
For an inhomogeneous atmosphere, the reflected light phase curve is still given by equation (12). The expression relating the secondary eclipse depth D to the geometric albedo A g is [7] Another way to understand the preceding expression is to realise that the flux of reflected starlight at secondary eclipse is 2F 0 I (the factor of 2 accounts for the fact that the expression for F assumes latitudinal symmetry), which is divided by the flux from a Lambertian disk πI to obtain A g .

Fit to data
We simultaneously fit the long-cadence fluxes of the hot Jupiter-host Kepler-7 observed by the Kepler mission, spanning the times corresponding to the phase curve and the secondary eclipse. We retrieve the Pre-search Data Conditioning Simple Aperture Photometry (PDCSAP) flux light curve from MAST with lightkurve [53], and mask 8σ outliers in flux. For the purposes of this analysis, we have masked out transits. There are 60,934 remaining fluxes at 30 minute cadence from Quarters 0-17. The light curve model consists of four components: reflected light, thermal emission, the secondary eclipse and a Gaussian process. All components are fit simultaneously to the entire out-oftransit light curve.
The reflected light component is modelled with an inhomogeneous atmosphere and the Henyey-Greenstein law of reflection. The planet has two hemispheres defined by the single-scattering albedos 0 ≤ ω 0 and ω = ω 0 + ω ≤ 1, with one global scattering asymmetry factor g that follows a Gaussian prior N (0, 0.25), where the darker surface is bounded by the local longitudes −π/2 < x 1 < x 2 < π/2. We note that the parameter pair (ω 0 , ω) shares similar constraints to the quadratic limb-darkening parameters [54] and thus adapt the Note: The index i refers to the components S, L and C.
same technique to the single-scattering albedos for efficient, uninformative sampling.
The thermal emission is computed by parameterising the temperature map with a generalised spherical harmonic basis [55]. The thermal emission model requires two free parameters: the spherical harmonic power in the temperature map in the m = = 1 mode C 11 , and f = f 0 (1 − A B ). The heat redistribution parameter f 0 is bounded between 1/4 (full redistribution) and 2/3 (no redistribution) [22]. The Bond albedo is given by A B = qA g , where A g is the geometric albedo integrated over the Kepler bandpass. Thus, A B is selfconsistently derived at each step in the chains from the other phase curve parameters. Since infrared phase curves have not been measured for Kepler-7 b, an empirical relationship between the peak offset of the thermal component and the irradiation temperature is used [56,57], which yields an expected peak offset consistent with zero; we fix the orbital phase of maximum thermal emission to α = 0.
The secondary eclipse is modelled with the exoplanet package with no limb-darkening [58,59]. The mid-transit time, period, impact parameter, exoplanetary radius and stellar radius [16], as well as the stellar density [60], are fixed to their previously reported values. The ratio of the exoplanetary flux to the stellar flux is assumed to drop to zero during secondary eclipse. We account for the 30 minute exposure times by super-sampling the eclipse curve.
Finally, we account for the stellar rotation and time-correlated systematics in the photometry with Gaussian process (GP) regression, using the simpleharmonic oscillator (SHO) kernel implemented in celerite2 [61,62]. We set the kernel undamped period and the damping timescale to stellar rotation period (15.3 days) and twice the rotation period, respectively, as measured by the autocorrelation function of the observations after transits and eclipses have been masked. This ensures that most of the periodic power in the GP is concentrated near periods similar to the stellar rotation or longer, with significantly less power near the 4.9 day orbital period of the planet.
The reflected light and thermal emission phase curve model components are implemented in theano [63] for compatibility with the exoplanet ecosys-  Metropolis-Hastings sampling is used to infer all of the previously mentioned parameters, as well as the uncertainties on the 30 minute cadence fluxes, and find a maximum likelihood uncertainty of 130 ppm. The posteriors are integrated until the Gelman-Rubin statistic becomesr 1.05 for all parameters [65]. Several other open-source software packages were used for this research [66,67,68,69,70,71,72]. We enumerate the maximum-likelihood parameters and other derived quantities in Table 2.