© The Authors 2024

1 Introduction

Optical phase curves of exoplanets can constrain properties of planetary atmospheres (see reviews by Deming & Seager 2017; Parmentier & Crossfield 2018). The phase curve encodes degenerate information about the temperature of the photosphere or surface, the reflectance or scattering properties of the photosphere or surface, and the effects of the planet on the host star such as ellipsoidal variations and Doppler beaming.

The amplitudes and shapes of reflected light phase curves are functions of fundamental properties of the photosphere or scattering surface (Cowan & Agol 2011; Demory et al. 2013; Hu et al. 2015; Shporer & Hu 2015; García Muñoz & Isaak 2015; Jansen & Kipping 2018; Farr et al. 2018; Mayorga et al. 2020; Luger et al. 2022; Fraine et al. 2021). These properties include the single-scattering albedo as a function of longitude, and the asymmetry of scattering events (forward or backward scattering). In general these quantities are also functions of wavelength, and broad-optical phase curves constrain the bandpass-integrated quantities. From these scattering parameters, the geometric albedo and integral phase function can be computed.

Reflected light phase curves of exoplanet atmospheres are contaminated by thermal emission, ellipsoidal variations, and Doppler beaming (e.g., Demory et al. 2011, 2013; Shporer et al. 2014; Esteves et al. 2015). The precise ratio of reflected light to thermal emission is often poorly constrained by optical observations alone, since the contribution from thermal emission often peaks at a similar orbital phase to the orbital phase of maximal reflected light for a homogeneous reflector (Demory et al. 2013). For short-period planets, the hottest point on the dayside of a planet with strong atmospheric circulation may not be at the substellar point, which can sometimes be equally described without circulation but with inhomogeneous cloud coverage, when measured with photometry in the optical and near-infrared. The ellipsoidal and Doppler contributions to the phase curve are constrained by the orbital period of the planet with uncertain amplitudes. Due to the uncertain mixture of these contributions in a given phase curve, all of these effects must be considered simultaneously.

The joint inference required for studying planets in reflected light is strongly motivated by the planetary constraints revealed by reflected light. Phase curves can directly constrain two optical properties of condensates in the atmosphere, namely the single scattering albedo and scattering asymmetry parameter, which are linked to the species of condensates responsible for the scattering (Kitzmann & Heng 2018). Furthermore, with the analytic reflection model in Heng et al. (2021), the shape of the reflected light phase curve can be linked to a scattering phase function. Common phase functions used in the exoplanet literature include Rayleigh (Hubbard et al. 2001), isotropic (Demory et al. 2013), Henyey-Greenstein (Sudarsky et al. 2000; de Kok & Stam 2012; Robinson 2017), and double Henyey-Greenstein (Heng & Li 2021; Adams et al. 2022) for example. The choice of a particular planetary scattering phase function can be directly linked to the scattering phase function of condensates, given a particle composition and size as well as cloud coverage. The Solar System planets offer guideposts with more complex scattering phase functions. Mayorga et al. (2016, see Fig. 5) and Dyudina et al. (2016, see Figs. 4, 5) found that flyby observations of Jupiter are inconsistent with pure isotropic or Rayleigh scattering phase functions¹. These precision Solar System observations and theoretical advances indicate that the scattering phase function can and should be inferred from the phase curves themselves.

Stellar rotation and instrumental systematics generate correlated signals in the aperture photometry of exoplanet host stars. A Bayesian retrieval framework for exoplanet phase curve parameters must also account for the uncertainty in the detrend-ing technique used to account for stellar rotation and instrumental systematics. Gaussian process (GP) regression is an effective and efficient technique for removing these time-correlated signals, and propagating the uncertainties in the kernel hyperparam-eters that describe the correlation in time into the atmospheric parameters retrieved for a system (see e.g. Angus et al. 2018).

Markov chain Monte Carlo (MCMC) methods are commonly used in astrophysics to estimate posterior distributions for uncertain parameters from measurements (see review by Sharma 2017). Recently, gradient-based inference techniques are becoming more widely used in photometric analyses of exoplanets, such as Hamiltonian Monte Carlo (HMC) and its popular variant, the No U-Turn Sampler (NUTS; see for example Luger et al. 2019; Agol et al. 2020; Colón et al. 2020; Price-Whelan et al. 2020; Stefánsson et al. 2020; Dalba et al. 2021; Daylan et al. 2021; Agol et al. 2021; Van Eylen et al. 2021; Foreman-Mackey et al. 2021). NUTS is especially useful for efficiently sampling high dimensional and degenerate posterior distributions (Betancourt 2017). Many of the probabilistic programming frameworks that implement HMC methods also provide methods for leave one out cross validation, which can be used to compare and select among a suite of models. This approach has been applied more recently in astronomical contexts (Welbanks et al. 2023; Gandhi et al. 2023; Challener et al. 2023).

In this work, we develop a Bayesian inference framework for determining reflected light properties of hot Jupiters observed with Kepler. We outline the photometric model and its sampler in Sect. 2. We present the results of the photometric analysis and a predictive accuracy/model comparison technique in Sect. 3, including assessments of inhomogeneity in exoplanet albedos, time variability, and inference of the single scattering phase function. We interpret these results in terms of potential condensates and observational biases in Sect. 4. We discuss the implications in Sect. 5, and conclude in Sect. 6.

2 Methods

In this section, we outline the numerical methods for sampling from the posterior distributions for the single-scattering albedos and scattering asymmetry parameters of exoplanet atmospheres with optical phase curve photometry. The source code is freely available in a package called kelp². For computational efficiency, the phase curve models have been implemented in JAX (Bradbury et al. 2018), a transpiled framework in Python. All phase curve parameter priors are listed in Table 1.

2.1 Photometry and system parameters

In this work, we focus our attention on planets with equilibrium temperatures T_eq ≳ 1500 K, as the daysides of these planets are likely to transition from cloudy to cloud-free as temperature increases above the stability limits of condensates near the sub-stellar longitude. We selected five planets spanning a range of equilibrium temperatures above 1500 K with geometries favorable for measuring reflected light. The sample of selected planets is shown in Fig. 1.

We fix the orbital period, mid-transit time, the stellar mass, planetary mass, planetary radius, semimajor axis, stellar density, impact parameter, and stellar temperature to literature values curated by the NASA Exoplanet Archive in the PSCompPars table.

We retrieve each long-cadence light curve from MAST as Simple Aperture Photomery (SAP) fluxes. We detrend the unmasked light curves with cotrending basis vectors (CBVs) also downloaded from MAST³ via lightkurve (Lightkurve Collaboration 2018), extracting the first eight basis vectors and applying single-scale CBV correction with the L2-norm regular-ization penalty α = 10⁻⁴.

Fig. 1 Host star Kepler magnitudes and their planetary equilibrium temperatures for the five planets we consider (large points). The color of each point indicates the amplitude of the reflected light phase curve assuming Ag = 1 (upper limit).

2.2 Reflected light

The reflected light phase curve is computed with the ab initio solutions for any reflection law by Heng et al. (2021). We test models with both homogeneous and inhomogeneous planetary atmospheres, where the latter is a generalization of the piecewise-Lambertian model of Hu et al. (2015). The inhomo-geneous sphere has a dark region with single scattering albedo 0 ≤ ω₀ ≤ 1 which is surrounded by a brighter region 0 < ω ≡ ω₀ + ω' < 1. The less reflective region is bounded by longitudes −π/2 < x₁ < x₂ < π/2. The geometric albedo is the last free parameter, which completes the set of parameters necessary to derive the scattering asymmetry parameter −1 < g < 1.

We apply a prior to the scattering asymmetry parameter g ~ 𝒩(0, 0.01) since the approximations in the reflected light model are most accurate for strongly asymmetric phase functions when g is closer to zero. The reflected light model for an inhomogeneous atmosphere is given by $$F_{\text{reflect }}={\left(\frac{R_{\text{p}}}{a}\right)}^2{A}_{\text{g}}\text{Ψ}\left({\omega }_0,\omega ,g,x_1,x_2\right)\text{, }$$(1)

where R_p is the planet radius, a is the semimajor axis, A_g is the geometric albedo, and the functional form of the integral phase function Ψ is stated in Eqs. (40) and (42) of Heng et al. (2021). When ω₀ = ω, the inhomogeneous sphere model reduces to a homogeneous sphere.

The less-reflective substellar region is bounded by longitudes x₁ and x₂, defined relative to the substellar longitude. To sample longitudes uniformly in the reference frame of a distant observer, we sample $x_1={\sin }^{-1}\left(x_1^{\prime }\right)\text{ and }{x}_2={\sin }^{-1}\left(x_2^{\prime }\right)$, with $x_1^{\prime }~\mathcal{U}(-1,0)\text{ and }{x}_2^{\prime }~\mathcal{U}(0,1)$.

2.2.1 Single scattering albedo

The two-albedo inhomogeneous model was inspired by the phase curve of Kepler-7 b (Hu et al. 2015). The physical motivation for the two albedos used for Kepler-7 b comes from theoretical expectations that cloudy regions have high albedos in comparison with low-albedo cloud-free regions. The temperature near the substellar longitude in these atmospheres is likely too hot for condensates to persist, though at some longitude closer to the limb, the temperatures may dip below the stability curves of highly refractory species. The inhomogeneous model captures this behavior with a dark substellar region with albedo ω₀, bounded by freely varying longitudes to the west and east where the albedo increases to ω.

Phase curve analyses with two-albedo models usually infer that the substellar regions have near-zero albedo (Hu et al. 2015; Adams et al. 2022). We confirm that fits to the phase curves with a uniform, uninformative prior on ω₀, and we find that all phase curves in this work are consistent with ω₀ ≈ 0. In all subsequent fits, including all results presented in this work, we simplify the inhomogeneous model parameterization by fixing ω₀ = 0 in the substellar region. In addition to improving fit convergence, fixing ω₀ simplifies the interpretation of the results in Sect. 4, since we can cleanly assume the region between x₁ and x₂ is cloud-free.

We note that the single-scattering albedo constrained in this work is defined mathematically by Heng et al. (2021), and its application to Kepler phase curves has some subtleties. The single-scattering albedo is defined for a particular wavelength. Since Kepler observations span a red-optical bandpass, the single-scattering albedos reported here are flux-weighted in the Kepler bandpass. The albedo measured here is also a weighted mean over the single-scattering albedo of the gas as well as the single-scattering albedo of the clouds.

Fig. 2 Scattering phase function proposed by Cornette & Shanks (CS, 1992), which reduces to Rayleigh when g = 0.

2.2.2 Scattering phase function

The phase curve formulae from Heng et al. (2021) allow users to trivially substitute different scattering phase functions for the exoplanet atmosphere. In principle, the choice of scattering phase function maps onto assumptions about the particle size in the planet’s atmosphere, the type of condensate (refractive indices), cloudiness, and the proportion of single-to-multiple scattering. Three simple scattering phase functions are considered in this work: (1) isotropic, which is convenient to implement, though practically one should not expect phase functions from single scattering events to be isotropic; (2) Rayleigh, which is expected for clear or cloudy atmospheres with particles that are much smaller than the wavelength of light considered (Strutt 1871, 1899); and (3) the Cornette & Shanks (1992) scattering phase function, $$P=\frac{3}{2}\frac{1-g^2}{2+g^2}\frac{1+{\mu }^2}{{\left(1+g^2-2g\mu \right)}^{3/2}},$$(2)

shown in Fig. 2. The latter is a generalization of the Rayleigh scattering phase function with additional freedom for stronger forward or backward scattering, parameterized by the scattering asymmetry parameter g. In the limit of g = 0, the Cornette-Shanks scattering phase function correctly reduces to Rayleigh instead of isotropic scattering.

2.3 Thermal emission

We model thermal emission of each exoplanet in the optical with the thermal emission model defined in Morris et al. (2022). The so-called h_mℓ basis describes the temperature map of the exo-planet using generalized spherical harmonics (parabolic cylinder functions, Heng & Workman 2014). There are two parameters which define the strength of the chevron shape of the hotspot feature on the dayside, which we fix to the values compatible with both GCM temperature maps and Spitzer phase curves, α = 0.6 and ω_drag = 4.5. Details and validation for this approach are provided in Morris et al. (2022).

We sample for two free parameters that define the atmosphere’s temperature map, which produces thermal emission. The first is the power in the m = ℓ = 1 spherical harmonic term, C₁₁, which determines the contrast between the dayside and nightside temperatures. Expected values for C₁₁ are of order 0.1, so we draw from a normal distribution C₁₁ ~ 𝒩(0, 0.1), with lower and upper limits at zero and 0.2. The lower limit corresponds to the limiting case of perfect heat redistribution from the dayside to the nightside, and the upper limit prevents extreme day/night contrasts.

The second free parameter is f, which is sometimes called the greenhouse factor. If we write the temperature map of the planet as a perturbation around a mean temperature $\overline{T}$, which is similar to the equilibrium temperature, $$\overline{T}=f{T}_{\star }\sqrt{R_{\star }/a},$$(3)

where f is related to the atmosphere’s efficiency at converting incoming stellar radiation into thermal emission (see Sects. 2 and 3 of Morris et al. 2022 for more discussion). We sample within a small range about the expected value f ~ 𝒩(2^−1/2, 0.01).

We find that the prior’s upper limit C₁₁ < 0.2 is necessary for fits to optical phase curves in this parameterization. Primarily, f sets the global mean temperature and C₁₁ sets the day/night contrast. Since the great majority of the flux from thermal emission is detected at full phase (ξ = 0), the hemisphere-averaged dayside temperature would be the best-constrained quantity to measure from an thermal phase curve in the optical. The dayside temperature in this parameterization results from a combination of f and C₁₁, where these terms are somewhat degenerate − a higher global mean temperature and smaller contrast can produce the same dayside temperature as a lower global mean temperature and a higher contrast. In infrared phase curves, as in Morris et al. (2022), this degeneracy can be avoided because the nightside thermal emission can be clearly detected. The nightside emission is negligible in the optical Kepler bandpass, so an upper limit on C₁₁ prevents sampling extreme solutions with hot daysides, cool nightsides, and low global mean temperatures.

We assume zero hotspot offset for each planet, ∆ϕ = 0. Spitzer phase curve analysis in Bell et al. (2021) and Morris et al. (2022, see their Fig. 8) showed that the hotspot offset is often small for planets with 1000 K < T_eq < 2500 K, which encompasses all planets in this sample. Practically, assuming zero hotspot offset also provides the maximum contamination of the reflected light signal from thermal emission, and allowing C₁₁ to vary allows us to explore the degeneracy between thermal emission and reflected light at secondary eclipse. Fixing ∆ϕ = 0 therefore implies that the geometric albedos that we report may be lower limits, since a non-zero hotspot offset will increase the contribution of reflected light to the phase curve amplitude.

2.4 Ellipsoidal variations and Doppler beaming

We use the following relations to estimate the amplitudes of the ellipsoidal variations and Doppler beaming respectively (see e.g. Shporer et al. 2014, and references therein): $$A_{\text{ellip }}[\text{ppm}]\approx \frac{{\alpha }_{\text{ellip }}}{0.077}\frac{M_{\text{p}}}{M_{\text{J}}}{\left(\frac{R_{\star }}{R_{\odot }}\right)}^3{\left(\frac{M_{\star }}{M_{\odot }}\right)}^{-2}{\left(\frac{P}{\text{ day }}\right)}^{-2}$$(4) $$A_{\text{Doppler }}[\text{ppm}]\approx \frac{{\alpha }_{\text{Doppler }}}{0.37}\frac{M_{\text{p}}}{M_{\text{J}}}{\left(\frac{M_{\star }}{M_{\odot }}\right)}^{-2/3}{\left(\frac{P}{\text{ day }}\right)}^{-1/3}$$(5)

where each α coefficient is a factor of order unity.

The flux from ellipsoidal variations has the form $$F_{\text{ellip }}=A_{\text{ellip }}\left(1-\cos 4\pi \left({\varphi }_{\text{orb }}-0.5\right)\right),$$(6)

where the orbital phase ϕ_orb is normalized to vary from zero at mid-transit to 0.5 at secondary eclipse.

The Doppler beaming model has the form $$F_{\text{Doppler }}=A_{\text{Doppler }}\sin 2\pi {\varphi }_{\text{orb }},$$(7)

and generally has a smaller amplitude than the ellipsoidal variations.

2.5 Eclipse model

The secondary eclipse model for each planet λ^e is the occupation defined by Agol et al. (2020) as implemented in exoplanet (Foreman-Mackey et al. 2021). We normalize λ^e such that it is unity out-of-eclipse and zero between second and third contacts, so that we can multiply the planetary contributions to the phase curve by this function.

2.6 Flux contamination from nearby stars

In Kepler observations, the flux from the target star may be diluted by flux from neighboring stars which falls within the photometric aperture of the target star. For each target except KOI-13, we adopt the Simple Aperture Photometry (SAP) crowding factor CROWDSAP in the Kepler and TESS input catalogs as the fraction of flux in the aperture from the target star. For KOI-13, we adopt the crowding factor due to flux from the other components in the hierarchical triple measured by Shporer et al. (2014).

2.7 Composite mean model

The complete phase curve model is given by $$\frac{F_{\text{p}}}{F_{\star }}=\left({\lambda }^e\left(F_{\text{reflect }}+F_{\text{thermal }}\right)+F_{\text{ellip }}+F_{\text{Doppler }}\right)/{\delta }_{\text{F}},$$(8)

where δ_F is the crowding factor as described in Sect. 2.6.

2.8 Gaussian process regression

We incorporate a Gaussian process regression step into the phase curve analysis using the Matérn 3/2 kernel implemented in celerite2 (Foreman-Mackey 2018). We use a fixed timescale ρ_GP = 30 days and unknown standard deviation of the process σ_Gp. The fixed ρ_GP timescale is about ten-times longer than the phase curve variations (P < 3 days), and thus removes instrumental and stellar rotation signals without affecting the phase curve.

2.9 Posterior sampling

We construct the composite mean model and Gaussian process marginal likelihood in the numpyro framework with JAX (Bradbury et al. 2018; Phan et al. 2019). The numpyro inference framework allows Python code to be compiled for efficient, parallel, Monte Carlo posterior sampling. numpyro provides support for several samplers, including standard Metropolis-Hastings and the No U-Turn Sampler (NUTS), which is a variation of the gradient-based Hamiltonian Monte Carlo (HMC) technique (for a review of HMC, see Betancourt 2017). We choose the No U-Turn Sampler for this work because in comparison with Metropolis-Hastings, for example, it often produces more effective samples from the posterior distribution in the less time⁴. For each exoplanet, we construct phase curve models with the homogeneous and inhomogeneous reflected light parameter-izations, and draw posterior samples with each model until the Gelman-Rubin statistic r̂ < 1.01 (Gelman & Rubin 1992; Vehtari et al. 2021).

2.10 Injection-recovery test

To evaluate the accuracy and precision of the kelp framework for retrieving planetary atmospheric properties, we inject a synthetic phase curve into the real Kepler photometry of KIC 8424992. This G2V solar analog is a reasonable target for phase curve injections, as it is similar to the mostly Sun-like stars in this sample, with K_p = 10.3 mag, without known planets, and the following asteroseismic stellar properties: mass 0.9 M_⊙, radius 1.05 R_⊙, and age 9.8 Gyr (Silva Aguirre et al. 2017).

To stress-test the framework, we inject a synthetic phase curve into the light curve of KIC 8424992 with significant albedo asymmetry, substantial thermal emission, and significant stellar artifacts. We inject an inhomogeneous reflected light phase curve similar in asymmetry to Kepler-7 b, and we give the planet thermal emission and orbital properties similar to HAT-P-7 b for significant thermal emission and stellar artifacts.

The posterior distributions for each of the stellar and planetary atmosphere parameters are shown in Fig. 3, along with the “true” (injected) values marked in blue. Most parameters are retrieved within 2σ of their input values. The least accurate measurements are the geometric albedo A_g (3.1σ), the start longitude of the darker region x₁ (2.8σ), and the spherical harmonic power C₁₁ (1.7σ), which are all degenerate with one other in the presence of ellipsoidal variations and Doppler beaming. These marginal inconsistencies arising from parameter degeneracies suggest that geometric albedos in this framework may have larger uncertainties than the posterior distributions suggest, up to a factor of a few.

3 Results

The maximum-likelihood models for all targets are shown in Fig. 4, and the reflected light and thermal emission fitting parameters are enumerated in Table 2. The four phase curve components are plotted both separately and combined. The silver points are the raw Kepler time series after applying the CBV and GP corrections. No binning in time was applied to these long-cadence (30 min) observations.

The time series after CBV correction but before phase curves and GPs are removed is shown in Fig. 5, in the time domain. Each light curve on the left shows 1 yr of Kepler observations, and their corresponding residual histogram is plotted in the matching color on the right. With the exceptions of infrequent outliers, the log-histogram on the right is approximately quadratic, which indicates nearly Gaussian noise.

3.1 Leave one out cross validation

After we estimate posterior distributions for the phase curve parameters, we use cross validation and stacking to select between the homogeneous and inhomogeneous reflected light models. Two information criteria often used are leave one out cross validation (LOO-CV) and the Widely Applicable (or Watanabe–Akaike) Information Criterion (WAIC). Following the recommendations in Vehtari et al. (2015a), we choose Pareto-Smoothed Importance Sampling (PSIS) LOO-CV for this analysis (Vehtari et al. 2015b). Since we construct our likelihood with a Gaussian process, we compute the LOO for a non-factorized model with the efficient algorithm of Bürkner et al. (2018). We briefly summarize this procedure below.

PSIS-LOO-CV can be constructed as an iterative procedure where the ith observation y_i is left out, denoted with the negative subscript y_−i to represent all observations except the ith. The log-likelihood of the remaining observations is recomputed, repeating this calculation for all observations. The sum of these log-likelihoods gives the expected log pointwise predictive density, elpd_LOO. When the likelihood we are evaluating is the marginal likelihood of a Gaussian process, the model cannot be trivially factorized, since the elpd_LOO will depend on the draws for the GP hyperparameters.

Bürkner et al. (2018) show that for invertible covariance matrix C, one can define two quantities, $$\begin{array}{c}{g}_i={\left[C^{-1}y\right]}_i,\\ {\overline{c}}_{ii}={\left[C^{-1}\right]}_{ii},\end{array}$$

where the ii subscripts denote diagonal elements, such that the log predictive density becomes $$\log p\left(y_i\mid y_{-i},\theta \right)=-\frac{1}{2}\log (2\pi )+\frac{1}{2}\log {\overline{c}}_{ii}-\frac{1}{2}\frac{g_i^2}{{\overline{c}}_{ii}}.$$(9)

Since this procedure involves inverting the full covariance matrix, it is expensive to compute, incurring computation times comparable to the integration for the posterior samples for these models. We perform this step “offline”, after posterior sampling is complete. To further accelerate LOO-CV computation times, we make another simplifying assumption. We assume that the log pointwise predictive density can be computed on local subsets of the observations and then concatenated to produce an approximate global lpd_LOO. This assumption is based on the finite autocorrelation timescale of the Matérn kernel. The observed fluxes are uncorrelated on timescales of several times ρ_GP, and therefore the influence of widely separated observations on the marginal likelihood should be negligible.

In practice, we choose each Kepler Quarter – three months of elapsed observing time – to be the local subset of observations on which to compute each lpd_LOO. The size of the covariance matrix that must be inverted in Eq. (9) would be the square of the full Kepler time series, of size ≈44 000² elements. Using the proposed subset LOO calculation reduces the size of the covariance matrix by a factor of 256 to a more computationally inexpensive ≈2700² elements. This smaller matrix can be inverted on a desktop computer without special memory requirements.

The batched, subset estimate of the LOO will be least accurate for observations at the beginning and end of Kepler quarters, on the ends of the batches. At these times, we will not propagate information into the marginal likelihood from autocorrelated astrophysical signals such as stellar rotation that are coherent across multiple Kepler quarters. Most stars considered in this work are solar-type stars and rarely have activity signals that are coherent for more than one month (Giles et al. 2017; Morris et al. 2019) – see Fig. 5. Thus, we posit that the subset LOO estimates may be sufficient for model selection in this case.

Fig. 3 Posterior correlation plot demonstrating an injection-recovery exercise with a synthetic inhomogeneous phase curve signal injected into the real Kepler light curve of KIC 8424992. The blue lines indicate the values chosen for each of the injected phase curve parameters. The most challenging inference in this analysis is determining which fraction of the phase curve signal comes from reflection or thermal emission, since higher geometric albedos and hotter daysides both increase the planetary flux. This produces the anti-correlation between C11 and Ag, since a cooler dayside can produce a similar phase curve to a more reflective atmosphere. The model produces accurate inferences despite this degeneracy by fitting for complementary parameters that define the shape of the phase curve at all phases, such as x1, x2, and g.

Fig. 4 phase curves of five hot Jupiters after CBV detrending and removing the Gaussian process, which captures stellar and instrumental artifacts. The four components of the phase curve are: reflected light (blue), thermal emission (orange), ellipsoidal variations (green), and Doppler beaming (dark red), along with the composite model in thick red. The unbinned 30-min Kepler time series is shown in gray circles and is used in the fit, and median-binned photometry is shown with black circles for ease of visualization.

3.2 Bayesian model selection

When several models are available to describe the observations, one may discard disfavored models by means of Bayesian model selection. Yao et al. (2018) recommend the stacking after comparing three types of Bayesian model averaging (BMA) with the model stacking technique. Both stacking and BMA methods yield values called the “model weights.” The model weights must sum to unity, and a model with weight close to zero can be rejected in favor of models with weights close to unity. Intermediate weights indicate no strong preference for either model.

We compute the model weights using a fork of the arviz.compare function. To confirm that stacking works within our framework, we verified that the PSIS-LOO-CV and stacking techniques indicate a clear preference for the inhomogeneous model when applied to the injection/recovery tests in Sect. 2.10, producing an inhomogeneous model weight close to unity for a phase curve with significant, prescribed inhomogeneity.

3.3 Homogeneity

Using the above LOO-CV and stacking techniques (see Sect. 3.1), we can efficiently compute the posterior predictive density for the full model including the GP, and determine which model best describes the observations. We compare models which assume the planet is an inhomogeneous reflector and models assuming homogeneous reflectors for each planet. The comparisons are plotted in Figs. 6 and 7 which shows darker blue points for planets more likely to be inhomogeneous reflectors. The results are also enumerated in Table 2.

The cross-validation only conclusively prefers the homogeneous phase curve model for TrES-2 b. The homogeneity selection is least confident about HAT-P-7 b and KOI-13 b (Table 2), since weights near 0.5 indicate no preference for either of the two models. We expect that the two least-certain homogeneity inferences occur for the two closest-in and hottest planets because they produce most thermal emission, and the largest ellipsoidal variations and Doppler beaming, all of which exacerbate degeneracies in the phase curve shape.

Theory also suggests less confident expectations for clouds on HAT-P-7 b and KOI-13 b. With dayside temperatures T_d > 2500 K, their atmospheres are likely too hot for condensates to persist, so we might expect them to be cloud-free, homogeneous reflectors (see for example Marley et al. 2013).

The inhomogeneity classification for HAT-P-7 b presents an interesting test case. With an equilibrium temperature of T_eq ≈ 2300 K, HAT-P-7 b resides close to the temperature above which we do not expect condensates to form. Dayside emission spectroscopy of HAT-P-7 b, like that of Mansfield et al. (2018), can be used to infer the temperature-pressure (p–T) structure for the atmosphere, which may be greater than this condensation limit (see for example Fig. 8 of Christiansen et al. 2010), but we emphasize here that these p-T profiles are typically a 1D dayside average, and do not account for the local surface temperature variations that we know must exist on real hot Jupiters. While the observations in this work indicate inhomogeneity, theoretical expectations for inhomogeneity are unclear this planet.

Kepler-7 b is perhaps one of the better known reflected light phase curves in the literature due to its high albedo, cool equilibrium temperature and strong reflected light asymmetry (Demory et al. 2013; Esteves et al. 2015; Hu et al. 2015; Angerhausen et al. 2015; Heng et al. 2021). As the target with the most clearly asymmetric phase curve, it is ideal for confirming that the model selection framework is performing correctly, and indeed the framework yields a strong preference for the inhomogeneous model. This may also be an important test of the GP flexibility – if the GP is too flexible, it could partially or entirely absorb the asymmetry in the phase curve of Kepler-7 b. It appears that the GP works as intended, and does not interfere with the short-period phase curve oscillations.

The confident detection of inhomogeneity for Kepler-41 b is in agreement with Shporer & Hu (2015), who pointed out that such an inference could be drawn from earlier works on the system (Santerne et al. 2011; Quintana et al. 2013). We comment on the results for Kepler-41 b in more detail in Sect. 4.1.

TrES-2 b is remarkable for its very low albedo (Kipping & Spiegel 2011). As the coolest planet in the sample, it should also have the least contamination from thermal emission in the Kepler bandpass. The 3σ upper limit on the geometric albedo is 2%. The LOO-CV analysis indicates a strong preference for inhomogeneity.

KOI-13 A is unique among the host stars for being a hierarchical triple, with the star the A component hosting the transiting planet (Shporer et al. 2014). We adopt the dilution due to the other components of the triple system from the analysis of Shporer et al. (2014), described in Sect. 2.6.

Fig. 5 Kepler SAP fluxes (gray points) after CBV corrections were applied, in relative flux units of ppm, for Quarter 2 (roughly 6% of each dataset). The composite fit (red curve) contains components including planetary phase curve (reflected light and thermal emission), strictly periodic stellar artifacts (ellipsoidal variations, Doppler beaming), and a Gaussian process for stochastic stellar artifacts such as magnetic activity (via a Gaussian process). The fit to the light curves on the left produces residuals which are shown in the histogram on the right.

3.4 Scattering phase functions

Heng et al. (2021) showed that the shape of the phase curve encodes information about the scattering phase function of the atmosphere. We examine here whether the phase curves contain sufficient information to identify the scattering phase function in each reflected light phase curve. Following Hapke (1981), Heng et al. (2021) assumes that multiple scattering is always isotropic, while the single-scattering component could be described by any phase function. We thus probe two independent questions in this subsection: (1) is isotropic multiple scattering a good assumption, and (2) do the observations prefer Rayleigh or Cornette & Shanks (1992) single scattering over isotropic.

Fig. 6 Inferred reflectance and temperature maps for each planet, where the center of the map corresponds to the substellar point. The shading in the upper row corresponds to the single scattering albedo ω of the reflective region on the dayside. Cross validation prefers inhomogeneous models for each planet, so each dayside features a dark sub-stellar region bounded by sub-observer longitudes x1 and x2· The thermal emission maps in the second row are not uniform, but appear nearly uniform in this figure because the colormap spans the full range of dayside temperatures for all five planets.

Fig. 7 Geometric albedos and equilibrium temperatures for the five Kepler planets, colored by the inhomogeneous model weight. Model comparisons with weights approaching unity confidently prefer the inhomogeneous model over the homogenous model (see Sect. 3.2), so darker points are more likely inhomogeneous reflectors. The marker shape corresponds to the preferred scattering phase function: squares are Cornette–Shanks single scattering plus isotropic multiple scattering, and circles are isotropic single plus isotropic multiple scattering. Results are also enumerated in Table 2.

3.4.1 Theoretical justification for isotropic multiple scattering

To test the applicability of the isotropic multiple scattering assumption made by Heng et al. (2021), we perform additional numerical simulations for Rayleigh and isotropic scattering phase functions. For the numerical calculations we use the discrete ordinate radiative transfer model C-DISORT (Hamre et al. 2013), which is the C-version of the well-established DISORT model (Stamnes et al. 1988).

DISORT provides the exact solution to the plane-parallel radiative transfer equation, taking into account thermal emission, surface scattering, and illumination by a beam source at the top of the atmosphere. It supports general scattering phase functions, described by the usual expansion into a Legendre series. Given a number of computational streams (ordinates), DISORT provides the full, angular-resolved intensity field I(µ, ϕ), where µ is the cosine of the polar and ϕ the azimuth angle.

In order to simulate a semi-infinite atmosphere (Chandrasekhar 1960) as done in the calculations presented in Heng et al. (2021), we set the optical depth at the bottom of the atmosphere to τ = 10⁶. The surface albedo is set to zero, but in principle plays no role for the resulting outgoing radiation due to the very high optical depth assumed at the lower boundary. Since we are only interested in the reflected light phase curve, we turn off the thermal emission term within DISORT and set the incident stellar beam as an exterior source.

DISORT calculates the radiation field in local coordinates rather than in the observer’s coordinate system as described in Heng et al. (2021). We therefore use Eq. (13) from Heng et al. (2021), which is originally from Sobolev (1975), to convert the observer-centric phase angle α, longitude Φ, and latitude Θ to the local coordinates for a given atmospheric column. This applies in particular to the local zenith angle µ_* of the stellar beam. From the DISORT calculations we then extract the intensity at the polar angle µ(Φ, Θ) and azimuth angle ϕ(Φ, Θ, α) scattered into the direction of the observer. The corresponding angles are again obtained by using Eq. (13) from Heng et al. (2021).

The radiative transfer is calculated for a range of values for longitude, latitude, and phase angle. The resulting intensities reflected toward the observer as a function of α, Φ, and Θ are transformed into the reflection coefficient ρ and then numerically integrated over Φ and Θ to yield the Sobolev fluxes F(α), following Eq. (5) from Heng et al. (2021), which is originally from Sobolev (1975). Finally, these fluxes are converted into the integral phase function $$\text{Ψ}(\alpha )=\frac{F(\alpha )}{F(\alpha =0)}.$$(10)

In Fig. 8, we compare the outcome for Ψ of numerical calculations performed with DisoRT with the analytical phase curve description from Heng et al. (2021) for Rayleigh and isotropic scattering phase functions. These calculations are identical to the scenario presented in Fig. 3a from Heng et al. (2021), assuming a constant single scattering albedo of ω = 0.5 throughout the atmosphere.

The analytic phase curves of Heng et al. (2021) are in excellent agreement with the numerical results. The assumption of isotropic multiple scattering is a good approximation for full-Rayleigh scattering. While the Rayleigh scattering phase function shows small deviations from an isotropic behavior, this is only important for single scattering events. When averaged over many multiple scattering events in a semi-infinite atmosphere, however, the phase function tends towards a more isotropic behavior.

Fig. 8 Comparing two simple scattering phase functions for a single-scattering albedo of ω = 0.5 with both single and multiple scattering (solid curves) and with different single scattering phase functions but forcing isotropic multiple scattering (dashed, white curves), as in Heng et al. (2021). The difference between "full" Rayleigh scattering and Rayleigh single plus isotropic multiple scattering is negligible.

3.4.2 Single scattering phase function

For all planets, using the LOO-CV and stacking techniques detailed above, the Cornette–Shanks or isotropic single scattering phase functions are clearly preferred, with Bayesian stacking model weights close to unity. The preference for isotropic single scattering is rather surprising at first, since no known atmospheric constituent, molecules, or potential condensates are known to produce a truly isotropic scattering phase function.

Molecules scatter radiation according to the Rayleigh scattering phase function. For condensates, assuming a spherical shape in the simplest case, Mie scattering is applicable (Mie 1908). When particles are small compared to the wavelength, the Mie scattering phase function simplifies to the Rayleigh one (Bohren & Huffman 1998). For example, a hypothetical homogeneously reflecting planet may be uniformly cloud-free, and reflected light from cloudless giant planets should be dominated by Rayleigh scattering from H₂ and He, which we expect to scatter with a Rayleigh phase function.

However, we should expect that multiple scattering is significant for several of these objects, as they have non-zero values of the single scattering albedo ω. As shown in the previous subsection, even pure Rayleigh single and Rayleigh multiple scattering will produce a phase curve similar to isotropic multiple scattering. The difference between Rayleigh and isotropic single scattering is maximal near quadrature but on the order of 10% in amplitude. We expect such subtle features in the shape of the phase curve to be degenerate with the other phase curve features such as the ellipsoidal variations in the presence of photometric noise. As a result, the preference for or against isotropic scattering could be seen as a weak test of the cloudiness of an atmosphere, though the clear asymmetry in the phase curves of a planet like Kepler-7 b is a much less ambiguous test of inhomogeneity.

In Heng et al. (2021), the Henyey-Greenstein (HG) scattering phase function was adopted for Kepler-7 b. For scattering asymmetry parameters g → 0, the HG function approaches an isotropic scattering phase function. The prior on g in Heng et al. (2021) therefore enforced that the scattering phase function was close to isotropic. This work uses a Cornette-Shanks phase function, so the strong prior g → 0 enforces a nearly Rayleigh phase function.

3.5 Variability

The best-fit phase curves are in Fig. 4 (thick red curves), which are comprised of a combination of reflected light, thermal emission, ellipsoidal variations, and Doppler beaming (thin curves). All models are fit to the unbinned, 30-min long cadence Kepler photometry with a Gaussian process in the time domain and CBV correction.

Figure 5 shows the CBV-corrected SAP fluxes from the kepler mission along with a sample model which has a Gaussian process with a mean model described by the four phase curve components: reflection, thermal emission, Doppler beaming and ellipsoidal variations. The light curves are separated into four segments and labeled with the year at the beginning of each year of observations. The residuals after the phase curve and Gaussian process model have been applied are plotted as a histogram on the right of Fig. 5, showing consistent, Gaussian normal distributions.

The large-amplitude signals in the CBV-corrected SAP fluxes in Fig. 5 may have stellar origins. Kepler-41, a G2V host, shows a consistent autocorrelated behavior in the flux residuals with amplitude 2 ppt at all times during the Kepler observations. GOV star TrES-2 shows possible rotational signals in the first and fourth year of Kepler observations. HAT-P-7 (F8V) and KOI-13 A (A5-7V) show 250 ppm variations which may be rotational or instrumental. In all cases, the Gaussian process accounts for this autocorrelated signal on timescales of tens of days, without significantly affecting phase curve signals on one-day timescales.

The geometric albedos we infer for each planet for each year are plotted in Fig. 9 and listed in Table 3. There are no significant departures from a constant geometric albedo for any of the four planets.

4 Interpretation

The techniques in this work advance our ability to characterize exoplanet atmospheres with photometry alone. In this section, we provide interpretation of the results presented in the previous section, assisted by condensation and Mie scattering calculations, to interpret what species might be responsible for inhomogeneous scattering. We also outline the observational biases introduced by the scattering parameters used in this work.

Fig. 9 Geometric albedos for four Kepler targets, measured in four 1-yr segments of the Kepler light curves.

4.1 Identifying potential scattering species

Theoretical expectations and the observational evidence outlined above point to inhomogeneous albedo distributions for all planets. For planets in this temperature range, we expect that differential condensation with longitude will produce inhomogeneous albedo distributions. At the cooler end the sample where T_eq < 2000 K, there are several refractory species that one might expect to condense out of regions in the dayside atmospheres (Marley et al. 2013).

Figure 10 shows stability curves for an atmosphere with [M/H] = 0.38. The curves indicate temperatures and pressures where condensation may occur, and are calculated with the chemical equilibrium model FASTCHEM⁵ (Stock et al. 2018, 2022; Kitzmann et al. 2023), using the JANAF tables (Chase 1998) to derive the mass action constants for selected condensates as described by Kitzmann et al. (2023). Some of the most stable high-temperature condensate species across all pressures are AI₂O₃ (corundum), CaTiO₃ (perovskite), MgAl₂O₄ (spinel), Mg₂SiO₄ (forsterite), MgSiO₃ (enstatite), and Fe(s) (see e.g. Fortney 2005 or Gail & Sedlmayr 2013). Each species condenses at temperatures T ≲ 2000 K at 1 bar, so they are all possible condensates for Kepler-7 b, Kepler-41 b and TreS-2 b. We can further narrow the list of possible condensates and their properties with two methods outlined below.

The first method is to predict the possible values of ω and g for each species using Mie theory, and compare them with the observations. Given the indices of refraction for a given species, the wavelength of observations, and the particle size, one can directly compute the expected value for ω and g using LX-MIE (Kitzmann & Heng 2018). The predictions are shown in Fig. 11 for five refractory species which may condense in these exoplanet atmospheres, as a function of particle size.

Figure 11 allows us to roughly constrain the sizes of condensate particles in the exoplanet atmospheres. All inhomogeneous planet have single scattering albedos ω ≳ 0.9. This provides a lower limit on the particle size since all compatible albedos occur for these scattering species when their particle radii r ≳ 0.01 µm. Similarly, we can place a rough upper limit on the particle size of condensates, since each phase curve can be reproduced with small g, excluding large particle sizes r ≳ 0.2 µm. The Mie calculations therefore weakly constrain the cloud particle radius between 10⁻²–10⁻¹ µm.

The Mie calculations disfavor AI₂O₃ (corundum) and Fe (iron) as the condensates responsible for reflection from Kepler-41 b and Kepler-7 b. These atmospheres have single scattering albedos ω ≳ 0.9, which is too large to be produced for AI₂O₃ or Fe particles of any radius.

Figure 11 is computed for spherical particles of a single size, and a distribution of particle sizes must be convolved with the ω curves for estimates of a nonuniform particle size distribution. The Mie calculations for ω are also upper limits on ω, since the true ω of a, for example, hydrogen-dominated atmosphere containing these species will be somewhat smaller due to gas absorption. Thus, Mie theory can exclude some species and weakly constrain the particle size.

Figure 12 shows another technique for identifying which condensates best match the inhomogeneous phase curves. The reflected light model has two free longitudinal parameters x\ and x₂, which indicate the longitudes that bound the low-albedo region of the planet. At longitudes x₁ < x < x₂ in the hot substellar region of the dayside, the atmosphere is too hot for condensates to form. These longitudes corresponding to the cloudless-to-cloudy transition are therefore also the longitudes where the local temperatures T(x₁) and T(x₂) drop below the stability curve of the cloud species. The thermal emission phase curve model, constructed with the h_mℓ basis functions, describes the temperature of the planet at each latitude and longitude, so we may calculate the local temperatures from the best-fit thermal emission map T(x₁) and T(x₂), evaluated at the longitudes of the transition from cloudless-to-cloudy regions x₁ and x₂ from the reflected light distribution. We can then compare the inferred transition longitudes and their corresponding temperatures to the stability curves of likely condensate species. If the inferred (x₁, x₂) coordinates from the reflected light model are consistent with a species’ predicted starting and ending longitudes, we have identified a candidate condensate for that atmosphere.

There is at least one more under-constrained parameter to consider in the calculation outlined above. Temperatures within planetary atmospheres vary with latitude and longitude, as we have encoded with the h_mℓ basis, but also with altitude or pressure, which we have neglected up to this point. Just as the condensation curves of various species are crossed as one moves farther from the sub-stellar point, temperature decreases with increasing altitude and therefore condensation can only occur above some altitudes. So theoretical predictions for the start/stop longitudes of the dark region are themselves functions of pressure, (x₁(p), x₂(p)). For each species, starting and ending longitudes of the dark region are thus compatible with some range of pressures where condensation may occur. The unknown cloud-deck pressure is not a fitting parameter in our models, and there are few theoretical expectations for a reasonable range, so we evaluate a range of pressures for each species.

The three panels of Fig. 12 show the joint posterior distribution for the starting and ending longitudes of the dark region, x₁ and x₂ in grayscale. The track of colored points joined by a silver line show the (x₁(p), x₂(p)) coordinates predicted from combining the stability curve of each species with the temperature map for Kepler-41 b. The posterior distribution for (x₁, x₂) is 2σ-consistent with the theoretical predictions for perovskite, forsterite, and enstatite clouds placed at pressures ranging from 0.08 to 0.2 bar.

In summary, theoretical considerations from Mie theory and condensation maps yield the following key inferences: (1) two of the five most refractory species can be ruled out by their scattering properties; (2) the plausible particle sizes for these condensates is 10^−2–10⁻¹ µm; (3) the three plausible remaining species could condense in the correct locations in the atmosphere to reproduce the inhomogeneity observed in the phase curves; (4) the pressure at the cloud layer may be in the range 0.08 to 0.2 bar.

Fig. 10 Stability curves of selected high-temperature condensates for an atmosphere with the metallicity of Kepler-41 b [M/H] = 0.38 (Bonomo et al. 2015). The calculations are assuming chemical equilibrium in the gas phase.

Fig. 11 Predictions from Mie theory for the reflectance properties of several refractory species assuming spherical, monodisperse particles observed in the Kepler bandpass as a function of particle radius (Kitzmann & Heng 2018, curves above), for comparison with the posterior distributions for the reflectance parameters from Kepler observations (histograms). The measured single scattering albedos ω of each planet (also shown in Fig. 14) have 0.7 ≲ ω ≲ 1, placing a lower limit on the particle radius r ≳ 0.01 µm. The phase curves are consistent with scattering asymmetry parameters g < 0.2, though posterior distributions are so narrow as a result of the prior distribution. Fits consistent with near-zero g correspond to upper limits on particle radii of r ≲ 0.1 µm. The three species consistent with the phase curve observations for Kepler-7 b and Kepler-41 b are CaTiO3 (perovskite), Mg2SiO4 (forsterite), and MgSiO3 (enstatite), given that they have ω > 0.9 and near-zero g. We compare the condensation temperatures of each of these species with the local temperatures at the longitudes of the less reflective region in Fig. 12.

Fig. 12 Joint posterior distributions for the planetary longitudes bounding the less reflective region (xı, x2) for Kepler-41 b in grayscale, in comparison with the theoretical predictions for several refractory species in each subplot, in color. The gray curves indicate the 1 and 2σ contours for (x1, x2). Each species’ condensation temperature is converted to a coordinate in (xı, x2), which varies as a function of the cloud pressure, producing the track of colored points which are labeled with the corresponding cloud pressure in units of bar. For example, the observed longitudes spanning the less reflective region of Kepler-41 b are 2σ consistent with CaTiO3 clouds located near 0.19 bar. Plausible cloud pressures for all species range from 0.08 to 0.2 bar.

Fig. 13 Relative contributions to the geometric albedo from Rayleigh single scattering events and isotropic multiple scattering (Heng et al. 2021).

4.2 Observational biases

4.2.1 Single scattering albedo

The first, trivial bias implied by parameterizing ω arises from the intrinsic brightness of the planet. A more reflective planet with larger ω will produce a larger signal in reflected light and can therefore be detected more easily.

A more subtle observational bias arises from the relative contributions of single and multiple scattering for different single scattering albedos. In Fig. 13, the curves represent the fraction of the geometric albedo contributed by Cornette-Shanks single scattering, which can be computed trivially with Eqs. (22) and (24) of Heng et al. (2021, see also their Fig. 2a). Single scattering accounts for more than 80% of the geometric albedo of the planet when ω < 0.4 for g ≈ 0, and the contribution from multiple scattering becomes non-negligible for ω > 0.4.

This suggests a slight limitation of the technique presented here and in Heng et al. (2021), since the framework can apply any single-scattering phase function to the reflected light phase curve, but the multiple scattering is assumed to be isotropic. Figure 13 shows that at large single scattering albedos, a significant fraction of the reflected light phase curve is contributed by isotropic multiple scattering. As a result, we should expect the shape of the phase curve to constrain the single scattering phase function for ω < 0.4; otherwise the assumption of isotropic multiple scattering begins to affect the shape of the reflected light phase curve. This partly explains why we found in Sect. 3.4.1 that the isotropic multiple scattering assumption has only a very small effect on the reflected light phase curve shape at ω = 0.5.

It is possible that for small enough ω, the single scattering will swamp the multiple scattering signal, but the first-mentioned observational bias against detecting faint objects will make it more difficult to measure the shapes of these phase curves. Therefore, we conclude that the isotropic multiple scattering assumption in the Heng et al. (2021) reflected light model likely has a small effect on this analysis.

These results are in good agreement with the Monte Carlo radiative transfer simulations of HD 189733 b by Lee et al. (2017), who found that Rayleigh-like scattering phase functions are necessary where single scattering dominates, and when multiple scattering dominates, isotropic phase functions are a good approximation.

4.2.2 Scattering asymmetry parameter

Another bias affects phase curve scattering parameter inference through the asymmetry parameter g. A value of g < 0 implies a net back-scattering atmosphere while g > 0 indicates predominantly forward scattering. A back-scattering atmosphere will reflect more light near secondary eclipse than an atmosphere, which is forward scattering, and therefore the reflected light phase curve amplitude will be larger for planets with back-scattering atmospheres. Forward scattering may be detectable by highly precise observations near transit (García Muñoz & Cabrera 2018).

To demonstrate the effect of this bias on retrieving phase curve scattering parameters, we can inject many reflected light phase curves into random normal noise with standard deviations of 50 ppm, with no thermal emission component from the planet. We solve for the best ω and g which fit the light curve with the L-BFGS-B minimizer, and plot the true and best-fit values for each parameter in Fig. B.1. The biases presented in this section are each reflected in the variance of the recovered parameters as a function of ω and g. As outlined in the previous subsection, the variance in ω is largest for small ω_true, since the phase curve amplitude scales with ω_true. As g becomes more negative, strong back-scattering amplifies the peak of the phase curve near secondary eclipse and the best-fit g becomes more precise. For g > 0, the forward-scattering pushes more flux into the planet’s atmosphere, reflecting less to the observer near secondary eclipse, reducing the precision on g. The precision on g ranges from < 1% for g < −0.1, about 3% at g = 0, and greater than 10% for g > 0.5.

The observational bias against confident detections of g > 0 affects our ability to detect reflected light from large particles. As shown in Fig. 11, small particles (r < 0.01 µm) tend to reflect with g → 0, regardless of species. The strong forward-scattering nature of larger particles, which greatly reduces the observed flux in the phase curve, implies that we are biased towards detecting reflected light phase curves from planets with clouds made of smaller particles. An interesting exception to this trend is iron, which slightly back-scatters light when the particle radius is near 0.1 µm.

5 Discussion

5.1 Correcting the single scattering albedo for Kepler-7 b

In comparing our results for Kepler-7 b with Heng et al. (2021), we found that we derive similar values for each fitting parameter except for the single scattering albedo of Kepler-7 b. We note here that we have discovered a typo in the code that produced the corner plot in Fig. 4b and the table in Extended Data Fig. 6 of Heng et al. (2021), which incorrectly reported the single scattering albedo for the more reflective region, ω′, and its sum with the albedo less reflective region ω = ω₀ + ω′. Heng et al. (2021) reported ω = 0.12 ± 0.05, whereas the corrected calculation is in fact much closer to unity, ω = 0.99 ± 0.01. The corrected values are in much better agreement with the analysis by Hu et al. (2015) for example, whose Table 1 shows an analogous quantity to ω₀ and ω′ written as r₀ and r₁. We list all literature values and our revisions for the albedo of Kepler-7 b in Table C.1.

5.2 Phase curve photometry as a cloud detector

We have presented evidence for clear and cloudy skies in exo-planet atmospheres based on time-series photometry alone in Sect. 3. Traditionally this inference has been done via spec-troscopy, either in emission or transmission. For targets where both observations are available, a combination of photometry and spectroscopy should be used to infer cloudiness, but by precisely testing the shape of a single-band, optical phase curve we can infer cloudiness.

The inhomogeneity inferred from a phase curve is a stronger test of the partial cloudiness of an atmosphere when compared with the effect of the scattering phase function on the shape of the phase curve, as outlined in Sect. 3.4, because we currently have no way of producing the inhomogeneity from a cloud-free atmosphere. Inhomogeneity produces asymmetries in the phase curve that are readily detected with high S/N, Kepler-like observations. Identifying the Rayleigh scattering expected for a clear H₂/He atmosphere from the phase curve shape alone is quite challenging since the effect of different scattering phase functions on the phase curve amplitude is small.

5.3 No cloudy-to-clear transition temperature

It is tempting to imagine that there is a transition in equilibrium temperatures above which atmospheres are too hot to form condensates and are always clear or cloud-free. In practice, there is likely a transition regime where condensates may form but they do or do not based on more complicated physics than is considered in this work. As more planets are discovered and their phase curve photometry fills in Fig. 7, we may find cloud-free, homogeneous atmospheres throughout this temperature range.

5.4 Condensate species and their sizes

In Sect. 4.1, we propose a particle size range consistent with the optical properties of the exoplanet atmospheres in reflected light. Lecavelier Des Etangs et al. (2008) find a similar particle radius range for enstatite grains in the atmosphere of HD 189733 b, and a similar size range was later invoked by Pont et al. (2013) for the same planet. Enstatite grains have also been invoked to fit interferometric spectroscopy of HD 206893 B (Kammerer et al. 2021). Wakeford & Sing (2015) considered transmission spectroscopy of cloudy planets with varying compositions and particle sizes, and found that the spectrum of HD 189733 b could be adequately fit by enstatite clouds with particle sizes ranging from 10⁻²–10⁻¹ µm, consistent with the interpretation of this work.

García Muñoz & Isaak (2015) produced synthetic phase curves for Kepler-7 b as a function of ω and g, assuming a double Henyey–Greenstein scattering phase function, and performed a χ² minimization to select the most-likely properties of the scattering condensates. The likely species and particle sizes enumerated in García Muñoz & Isaak (2015) are in good agreement with our results from the case study on Kepler-41 b.

Webber et al. (2015) studied reflected light phase curves with inhomogeneous cloud cover applied to Kepler-7 b. The authors find that enstatite or forsterite condensation provides the best match to the observed reflected light asymmetry and albedo for Kepler-7 b. Similarly, Oreshenko et al. (2016) suggest that enstatite or forsterite (plus iron, corundum, or titanium oxide) provide reasonable matches to the phase offset observed for Kepler-7 b. The analysis presented in Sect. 4.1 can rule out some of the species considered in these previous works based on the scattering properties (ω, g) of the condensates.

Parmentier et al. (2016) used GCMs, condensation temperatures, and scattering properties of candidate condensates to constrain ensemble cloud properties of hot Jupiters. In the equilibrium temperature range of planets discussed in this work, the authors highlight silicate and perovskite clouds for cooler atmospheres (Kepler-7 b and Kepler-41 b, respectively) and corundum clouds for hotter atmospheres (HAT-P-7 b) as potential scatter-ers consistent with phase curve observations. Their simulated planets with the highest albedos and largest reflected light asymmetries are produced for particles with sizes near 10⁻¹ µm, which is consistent with the size upper limits implied by the small scattering asymmetries in this work.

Our results are also consistent with Lee et al. (2016, 2017), who presented 3D radiative-hydrodynamic simulations with a kinetic, microphysical, nonequilibrium mineral cloud model for HD 189733 b. Their condensate particle sizes are typically sub-micron, with cloud formation at the western limb and at mid-latitudes dominated by enstatite and forsterite.

Cloud microphysics models by Powell et al. (2019) suggest that the west limbs of hot Jupiters with T_eq = 2000 K may have, for example, forsterite condensation with number densities which peak near 10⁻² and 1 µm (see their Fig. 5). These particle sizes are consistent with the observations presented here, in the interpretation in Sect. 4.1. More generally, Gao et al. (2020) provide observational and theoretical evidence that silicates are the dominant condensates in hot atmospheres with T_eq ≳ 1600 K, which includes all planets in this sample.

5.5 Variability in exoplanet atmospheres

Uncertainty in time-series detrending should inform posterior inferences. The impact of the detrending techniques on the inferred phase curve parameters is a key challenge in interpreting reflected light phase curves in the literature. Often, conservative detrending approaches are applied independently from the posterior inference, producing, for example, phase curves that appear highly variable with time (Armstrong et al. 2016); while more liberal detrending risks removing signals of interest. Theoretically, there are few or no predictions for large variations in the thermal emission components of the phase curve on year-long timescales. Komacek & Showman (2020) used general circulation models to show that the variations in eclipse depth, phase curve amplitude, and hotspot offset should all change by a few percent or less. High resolution pseudospectral simulations by Y-K. Cho et al. (2021) are also in agreement with this amplitude of variations in the thermal emission. Larger variations might be expected in reflected light since, for example, there may be significant albedo variations driven by partial cloud cover, as we see on planets in the Solar System. The observations presented in this work do not support larger than a few percent in the geometric albedo of hot Jupiters on year-long timescales.

Claims of variations of exoplanet atmospheric properties should be approached with caution. Many observations support the theoretical expectations against significant thermal variability: Agol et al. (2010); Kilpatrick et al. (2020) found that the thermal flux of HD 189733 b and HD 209458 b change by less than a few percent between eclipses. Jones et al. (2022) found a similar upper limit on the optical flux variations of the ultra-hot Jupiter KELT-9 b. Recent work by Lally & Vanderburg (2022) showed that apparent variations like those invoked for the Kepler phase curve of HAT-P-7 b may be explained by uncorrected stellar or instrumental artifacts.

Fig. 14 Posterior distributions for the single scattering albedos. For inhomogeneous atmospheres, ω is the single scattering albedo of the more reflective region towards the limbs, and for homogeneous atmospheres ω is the dayside hemisphere-averaged single scattering albedo. For comparison, the single scattering albedos of plausible condensates are shown in Fig. 11.

6 Conclusion

• We have introduced a Bayesian inference framework for inferring optical properties of exoplanet atmospheres from phase curves, which includes flux contributions from reflected light from a potentially inhomogeneous atmosphere, thermal emission, ellipsoidal variations, Doppler beaming, and stellar rotation (Sect. 2);

• We applied this inference framework to the Kepler phase curves of five hot Jupiters to measure atmospheric homogeneity and time-variability, as well as an investigation into the scattering properties which constrain the likely condensates in inhomogeneous atmospheres. We use the leave one out cross validation technique with Bayesian stacking for model selection (Sects. 3.1–3.2);

• The observations suggest an inhomogeneous albedo distribution for three of the five planets (Fig. 7), which we interpret as asymmetric cloud distributions;

• We enumerate the phase curve parameters of each planet (Table 2), and show the single scattering albedos for the reflective region of each planet (Fig. 14);

• None of the planets exhibit significant geometric albedo variations in time (Fig. 9);

• For Kepler-41 b, we have identified perovskite, forsterite, and enstatite as possible scattering species consistent with the reflected light phase curves using condensate stability curves (Fig. 10) and Mie theory (Fig. 11). The condensate particle radii may be in the range 10⁻²–10⁻¹ µm, and the pressure at the clouds may be in the range 0.08–0.2 bar (Fig. 12);

• We demonstrated that analytic phase curves with isotropic multiple scattering are in excellent agreement with full Rayleigh multiple scattering calculations (Fig. 8).

Acknowledgements

We are grateful to our referee, Nicolas Cowan. We gratefully acknowledge the open source software which made this work possible: astropy (Astropy Collaboration 2013, 2018, 2022), ipython (Pérez & Granger 2007), numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), matplotlib (Hunter 2007), JAX (Bradbury et al. 2018), arviz (Kumar et al. 2019), numpyro (Phan et al. 2019), FASTCHEM (Stock et al. 2018, 2022; Kitzmann et al. 2023), LX-MIE (Kitzmann & Heng 2018), celerite2 (Foreman-Mackey et al. 2017; Foreman-Mackey 2018) exoplanet (Foreman-Mackey et al. 2021), lightkurve (Lightkurve Collaboration 2018), corner (Foreman-Mackey 2016), kelp (Morris et al. 2022). This research has made use of the SVO Filter Profile Service (http://svo2.cab.inta-csic.es/theory/fps/) supported from the Spanish MINECO through grant AYA2017-84089. This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program.

Appendix A Index of symbols

Table A.1 describes the symbols used in this work and their definitions.

Appendix B Biases

An illustrative battery of phase curve injection/recovery tests are shown in Figure B.1. These exercises demonstrate the biases which affect retrieval of the single scattering albedo ω and scattering asymmetry parameter g.

Fig. B.1 Comparison of the true (injected) single scattering albedo ω and scattering asymmetry parameter g with their best-fit values in simulated light curves of reflected light phase curves, with no thermal emission and Gaussian noise. As ωtrue → 0, the reflected light signal diminishes, reducing the precision in ωfit. For g > 0, the planet’s atmosphere scatters more light into the atmosphere than back at the observer, reducing the reflected light signal, and reducing the precision on gfit. The color of each point in the g panel denotes the value of ω for that point, and vice versa. Note that g is least precise for the smallest ω values (purple).

Appendix C Results for Kepler-7 b

There are many results for the geometric albedo in the literature, which we collate in Table C.1.

References

All Tables

All Figures

	Fig. 1 Host star Kepler magnitudes and their planetary equilibrium temperatures for the five planets we consider (large points). The color of each point indicates the amplitude of the reflected light phase curve assuming Ag = 1 (upper limit).
In the text

	Fig. 2 Scattering phase function proposed by Cornette & Shanks (CS, 1992), which reduces to Rayleigh when g = 0.
In the text

Fig. 3 Posterior correlation plot demonstrating an injection-recovery exercise with a synthetic inhomogeneous phase curve signal injected into the real Kepler light curve of KIC 8424992. The blue lines indicate the values chosen for each of the injected phase curve parameters. The most challenging inference in this analysis is determining which fraction of the phase curve signal comes from reflection or thermal emission, since higher geometric albedos and hotter daysides both increase the planetary flux. This produces the anti-correlation between C11 and Ag, since a cooler dayside can produce a similar phase curve to a more reflective atmosphere. The model produces accurate inferences despite this degeneracy by fitting for complementary parameters that define the shape of the phase curve at all phases, such as x1, x2, and g.

In the text

Fig. 4 phase curves of five hot Jupiters after CBV detrending and removing the Gaussian process, which captures stellar and instrumental artifacts. The four components of the phase curve are: reflected light (blue), thermal emission (orange), ellipsoidal variations (green), and Doppler beaming (dark red), along with the composite model in thick red. The unbinned 30-min Kepler time series is shown in gray circles and is used in the fit, and median-binned photometry is shown with black circles for ease of visualization.

In the text

Fig. 5 Kepler SAP fluxes (gray points) after CBV corrections were applied, in relative flux units of ppm, for Quarter 2 (roughly 6% of each dataset). The composite fit (red curve) contains components including planetary phase curve (reflected light and thermal emission), strictly periodic stellar artifacts (ellipsoidal variations, Doppler beaming), and a Gaussian process for stochastic stellar artifacts such as magnetic activity (via a Gaussian process). The fit to the light curves on the left produces residuals which are shown in the histogram on the right.

In the text

Fig. 6 Inferred reflectance and temperature maps for each planet, where the center of the map corresponds to the substellar point. The shading in the upper row corresponds to the single scattering albedo ω of the reflective region on the dayside. Cross validation prefers inhomogeneous models for each planet, so each dayside features a dark sub-stellar region bounded by sub-observer longitudes x1 and x2· The thermal emission maps in the second row are not uniform, but appear nearly uniform in this figure because the colormap spans the full range of dayside temperatures for all five planets.

In the text

Fig. 7 Geometric albedos and equilibrium temperatures for the five Kepler planets, colored by the inhomogeneous model weight. Model comparisons with weights approaching unity confidently prefer the inhomogeneous model over the homogenous model (see Sect. 3.2), so darker points are more likely inhomogeneous reflectors. The marker shape corresponds to the preferred scattering phase function: squares are Cornette–Shanks single scattering plus isotropic multiple scattering, and circles are isotropic single plus isotropic multiple scattering. Results are also enumerated in Table 2.

In the text

Fig. 8 Comparing two simple scattering phase functions for a single-scattering albedo of ω = 0.5 with both single and multiple scattering (solid curves) and with different single scattering phase functions but forcing isotropic multiple scattering (dashed, white curves), as in Heng et al. (2021). The difference between "full" Rayleigh scattering and Rayleigh single plus isotropic multiple scattering is negligible.

In the text

	Fig. 9 Geometric albedos for four Kepler targets, measured in four 1-yr segments of the Kepler light curves.
In the text

	Fig. 10 Stability curves of selected high-temperature condensates for an atmosphere with the metallicity of Kepler-41 b [M/H] = 0.38 (Bonomo et al. 2015). The calculations are assuming chemical equilibrium in the gas phase.
In the text

Fig. 11 Predictions from Mie theory for the reflectance properties of several refractory species assuming spherical, monodisperse particles observed in the Kepler bandpass as a function of particle radius (Kitzmann & Heng 2018, curves above), for comparison with the posterior distributions for the reflectance parameters from Kepler observations (histograms). The measured single scattering albedos ω of each planet (also shown in Fig. 14) have 0.7 ≲ ω ≲ 1, placing a lower limit on the particle radius r ≳ 0.01 µm. The phase curves are consistent with scattering asymmetry parameters g < 0.2, though posterior distributions are so narrow as a result of the prior distribution. Fits consistent with near-zero g correspond to upper limits on particle radii of r ≲ 0.1 µm. The three species consistent with the phase curve observations for Kepler-7 b and Kepler-41 b are CaTiO3 (perovskite), Mg2SiO4 (forsterite), and MgSiO3 (enstatite), given that they have ω > 0.9 and near-zero g. We compare the condensation temperatures of each of these species with the local temperatures at the longitudes of the less reflective region in Fig. 12.

In the text

Fig. 12 Joint posterior distributions for the planetary longitudes bounding the less reflective region (xı, x2) for Kepler-41 b in grayscale, in comparison with the theoretical predictions for several refractory species in each subplot, in color. The gray curves indicate the 1 and 2σ contours for (x1, x2). Each species’ condensation temperature is converted to a coordinate in (xı, x2), which varies as a function of the cloud pressure, producing the track of colored points which are labeled with the corresponding cloud pressure in units of bar. For example, the observed longitudes spanning the less reflective region of Kepler-41 b are 2σ consistent with CaTiO3 clouds located near 0.19 bar. Plausible cloud pressures for all species range from 0.08 to 0.2 bar.

In the text

	Fig. 13 Relative contributions to the geometric albedo from Rayleigh single scattering events and isotropic multiple scattering (Heng et al. 2021).
In the text

	Fig. 14 Posterior distributions for the single scattering albedos. For inhomogeneous atmospheres, ω is the single scattering albedo of the more reflective region towards the limbs, and for homogeneous atmospheres ω is the dayside hemisphere-averaged single scattering albedo. For comparison, the single scattering albedos of plausible condensates are shown in Fig. 11.
In the text

Fig. B.1 Comparison of the true (injected) single scattering albedo ω and scattering asymmetry parameter g with their best-fit values in simulated light curves of reflected light phase curves, with no thermal emission and Gaussian noise. As ωtrue → 0, the reflected light signal diminishes, reducing the precision in ωfit. For g > 0, the planet’s atmosphere scatters more light into the atmosphere than back at the observer, reducing the reflected light signal, and reducing the precision on gfit. The color of each point in the g panel denotes the value of ω for that point, and vice versa. Note that g is least precise for the smallest ω values (purple).

In the text