Observations of scattered light from exoplanet atmospheres

Optical phase curves of hot Jupiters can reveal global scattering properties. We implement a Bayesian inference framework for optical phase curves with flux contributions from: reflected light from a potentially inhomogeneous atmosphere, thermal emission, ellipsoidal variations, Doppler beaming, and stellar rotation via a Gaussian process in the time domain. We probe for atmospheric homogeneity and time-variability using the reflected light inferences for highly precise Kepler light curves of five hot Jupiters. We also investigate the scattering properties which constrain the most likely condensates in the inhomogeneous atmospheres. Cross validation prefers inhomogeneous albedo distributions for Kepler-7 b and Kepler-41 b, and a weak preference for inhomogeneity for KOI-13 b. None of the five planets exhibit significant variations in geometric albedo on one-year timescales, in agreement with theoretical expectations. We show that analytic reflected light phase curves with isotropic multiple scattering are in excellent agreement with full Rayleigh multiple scattering calculations, allowing for accelerated and analytic inference. In a case study of Kepler-41 b, we identify perovskite, forsterite, and enstatite as possible scattering species consistent with the reflected light phase curves, with condensate particle radii in the range 0.01-0.1 micron.


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. 2021;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 bandpassintegrated 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. 2011Demory et al. , 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 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 2018a).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 Figure 5) and Dyudina et al. (2016, see Figures 4, 5) found that flyby observations of Jupiter are inconsistent with pure isotropic or Rayleigh scattering phase functions1 .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 detrending 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 hyperparameters 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. 2021a).NUTS is especially useful for efficiently sampling high dimensional and degenerate posterior distributions (Betancourt 2017).Many of the probabilistic programming frameworks which 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 Section 2. We present the results of the photometric analysis and a predictive accuracy/model comparison technique in Section 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 Section 4. We discuss the implications in Section 5, and conclude in Section 6.

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

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

Parameter
Prior dist./valueInterval Reflected light  increases above the stability limits of condensates near the substellar 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 Figure 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 MAST3 via lightkurve (Lightkurve Collaboration et al. 2018), extracting the first eight basis vectors and applying single-scale CBV correction with the L2-norm regularization penalty α = 10 −4 .

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 inhomogeneous sphere has a dark region with single scattering albedo 0 ≤ ω 0 ≤ 1 which is surrounded by a brighter region 0 < ω ≡ ω 0 + ω ′ < 1.The less reflective region is bounded by longitudes −π/2 < x 1 < x 2 < π/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 ∼ N(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 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 equations ( 40) and (42) of Heng et al. (2021).When ω 0 = ω, the inhomogeneous sphere model reduces to a homogeneous sphere.The less-reflective substellar region is bounded by longitudes x 1 and x 2 , defined relative to the substellar longitude.To sample longitudes uniformly in the reference frame of a distant observer, we sample x 1 = sin −1 (x ′ 1 ) and

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 ω 0 , 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 ω 0 , and we find that all phase curves in this work are consistent with ω 0 ≈ 0. In all subsequent fits, including all results presented in this work, we simplify the inhomogeneous model parameterization by fixing ω 0 = 0 in the substellar region.In addition to improving fit convergence, fixing ω 0 simplifies the interpretation of the results in Section 4, since we can cleanly assume the region between x 1 and x 2 is cloudfree.
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 singlescattering albedo is defined for a particular wavelength.Since Kepler observations span a red-optical bandpass, the singlescattering 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 singlescattering albedo of the clouds.

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(Strutt , 1899)); and (3) the Cornette & Shanks (1992) scattering phase function, shown in Figure 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.

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 exoplanet 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 11 , which determines the contrast between the dayside and nightside temperatures.Expected values for C 11 are of order 0.1, so we draw from a normal distribution C 11 ∼ N(0, 0.1), Article number, page 3 of 22 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 T , which is similar to the equilibrium temperature, f is related to the atmosphere's efficiency at converting incoming stellar radiation into thermal emission (see Sections 2 and 3 of Morris et al. 2022 for more discussion).We sample within a small range about the expected value f ∼ N(2 −1/2 , 0.01).We find that the prior's upper limit C 11 < 0.2 is necessary for fits to optical phase curves in this parameterization.Primarily, f sets the global mean temperature and C 11 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 11 , 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 11 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 Figure 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 11 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.

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 Doppler [ppm] ≈ α Doppler 0.37 where each α coefficient is a factor of order unity.
The flux from ellipsoidal variations has the form 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 and generally has a smaller amplitude than the ellipsoidal variations.

Eclipse model
The secondary eclipse model for each planet λ e is the occultation defined by Agol et al. (2020) as implemented in exoplanet (Foreman-Mackey et al. 2021a).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.

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

Composite mean model
The complete phase curve model is given by where δ F is the crowding factor as described in Section 2.6.

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 d), and thus removes instrumental and stellar rotation signals without affecting the phase curve.For each exoplanet, we construct phase curve models with the homogeneous and inhomogeneous reflected light parameterizations, and draw posterior samples with each model until the Gelman-Rubin statistic r < 1.01 (Gelman & Rubin 1992;Vehtari et al. 2019).

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 Figure 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 1 (2.8σ), and the spherical harmonic power C 11 (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.

Results
The maximum-likelihood models for all targets are shown in Figure 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 minute) observations.
The time series after CBV correction but before phase curves and GPs are removed is shown in Figure 5, in the time domain.Each light curve on the left shows one year 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.

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 i th observation y i is left out, denoted with the negative subscript y −i to represent all observations except the i th .The loglikelihood of the remaining observations is recomputed, repeating this calculation for all observations.The sum of these loglikelihoods 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 can not 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, where the ii subscripts denote diagonal elements, such that the log predictive density becomes 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 Equation 9 would be the square of the full Kepler time series, of size ≈ 44,000 2 elements.
Article number, page 5 of 22 A&A proofs: manuscript no.main  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 C 11 and A g , 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 x 1 , x 2 , and g.
Using the proposed subset LOO calculation reduces the size of the covariance matrix by a factor of 256 to a more computationally inexpensive ≈ 2,700 2 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 quar-ters, 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;    We compute the model weights using a fork of the arviz.comparefunction.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 Section 2.10, producing an inhomogeneous model weight close to unity for a phase curve with significant, prescribed inhomogeneity.

Homogeneity
Using the above LOO-CV and stacking techniques (see Sect.Table 2. Reflected light and thermal emission phase curve parameters inferred from Kepler observations.For prior distributions, see Table 1.The first section gives results from model comparison with leave-one-out cross-validation (LOO-CV) and model stacking.Models are marked "Hom." when LOO-CV prefers the homogeneous reflected light model, "Inhom."when LOO-CV prefers the inhomogeneous model.The scattering phase function is included when preferred by LOO-CV with the same criteria.The posterior distributions for the reflected light and thermal emission parameters are in the second section, excluding g which was consistent with the prior for all planets.The derived quantities include the phase offset of the reflected light phase curve, the amplitudes of the reflected light and thermal emission components, and the dayside and nightside temperatures.See Section 2.10 for discussion on the uncertainties derived from these degenerate posterior distribution.

Homogeneity
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 comparison is plotted in Figure 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 ellip-soidal 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 Article number, page 9 of 22   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 Article number, page 10 of 22 Brett M. Morris et al.: Observations of scattered light from exoplanet atmospheres 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 shortperiod 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 Section 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 Section 2.6.

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?

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

Rayleigh Isotropic
Fig. 9. Comparing two simple scattering phase functions for a singlescattering 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.
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 Equation ( 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 Equation ( 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 towards 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 Equation (5) from Heng et al. (2021), which is originally from Sobolev (1975).Finally, these fluxes are converted into the integral phase function In Figure 9, 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 Figure 3(a) 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.

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

Variability
The best-fit phase curves are in Figure 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-minute 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 Figure 5, showing consistent, Gaussian normal distributions.
The large-amplitude signals in the CBV-corrected SAP fluxes in Figure 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.G0V 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 oneday timescales.
The geometric albedos we infer for each planet for each year are plotted in Figure 10 and listed in Table 3.There are no significant departures from a constant geometric albedo for any of the four planets.
lations, 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.

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 8 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 FastChem5 (Stock et al. 2018(Stock et al. , 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 hightemperature condensate species across all pressures are Al 2 O 3 (corundum), CaTiO 3 (perovskite), MgAl 2 O 4 (spinel), Mg 2 SiO 4 (forsterite), MgSiO 3 (enstatite), and Fe(s) (see e.g.Fortney 2005or 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 2018b).The predictions are shown in Figure 12 for five refractory species which may condense in these exoplanet atmospheres, as a function of particle size.
Figure 12 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 −2 − 10 −1 µm.
The Mie calculations disfavor Al 2 O 3 (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 Al 2 O 3 or Fe particles of any radius.
Figure 12 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 non-uniform particle size distribution.The Mie calculations for ω are also upper limits on ω, since the true ω of a, e.g., 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 13 shows another technique for identifying which condensates best match the inhomogeneous phase curves.The reflected light model has two free longitudinal parameters x 1 and x 2 , which indicate the longitudes that bound the low-albedo region of the planet.At longitudes x 1 < x < x 2 in the hot substellar region of the dayside, the atmosphere is too hot for condensates to form.These longitudes corresponding to the cloudless-tocloudy transition are therefore also the longitudes where the local temperatures T (x 1 ) and T (x 2 ) 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 1 ) and T (x 2 ), evaluated at the longitudes of the transition from cloudless-to-cloudy regions x 1 and x 2 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 1 , x 2 ) 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 Article number, page 13 of 22  11) 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 CaTiO 3 (perovskite), Mg 2 SiO 4 (forsterite), and MgSiO 3 (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 Figure 13.longitudes of the dark region are themselves functions of pressure, (x 1 (p), x 2 (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 clouddeck 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 Figure 13 show the joint posterior distribution for the starting and ending longitudes of the dark region, x 1 and x 2 in grayscale.The track of colored points joined by a silver line show the (x 1 (p), x 2 (p)) coordinates predicted from combining the stability curve of each species with the temperature map for Kepler-41 b.The posterior distribution for (x 1 , x 2 ) 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 −1 µ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.

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 Figure 14, the curves represent the fraction of the geometric albedo contributed by Cornette-Shanks single scattering, which can be computed trivially with Equations 22 and 24 of Heng et al. (2021, see also their Figure 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 14 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 Section 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.

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 backscattering 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 Figure 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 Figure 12, small particles (r < 0.01 µm) tend to reflect with g → 0, regardless of species.The strong forwardscattering 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.

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 Figure 4b and the table in Extended Data Figure 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 ω ≡ ω 0 + ω ′ .Heng et al. (2021) reported ω = 0.12 ± 0.05, whereas the corrected calculation is Article number, page 15 of 22 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 ω 0 and ω ′ written as r 0 and r 1 .We list all literature values and our revisions for the albedo of Kepler-7 b in Table C.1.

Phase curve photometry as a cloud detector
We have presented evidence for clear and cloudy skies in exoplanet atmospheres based on time-series photometry alone in Section 3. Traditionally this inference has been done via spectroscopy, 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 Section 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 2 /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.

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 Figure 7, we may find cloud-free, homogeneous atmospheres throughout this temperature range.

Condensate species and their sizes
In Section 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 −2 − 10 −1 µ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 χ 2 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 Section 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 scatterers 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 −1 µ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. (2016Lee et al. ( , 2017)), who presented 3D radiative-hydrodynamic simulations with a kinetic, microphysical, non-equilibrium mineral cloud model for HD 189733 b.Their condensate particle sizes are typically submicron, with cloud formation at the western limb and at midlatitudes 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 −2 and 1 µm (see their Figure 5).These particle sizes are consistent with the observations presented here, in the interpretation in Section 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.

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

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 (Section 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-oneout cross-validation technique with Bayesian stacking for model selection (Section 3.1-3.2). -The observations suggest an inhomogeneous albedo distribution for three of the five planets (Figure 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 (Figure 11) -None of the planets exhibit significant geometric albedo variations in time (Figure 10).-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 (Figure 8) and Mie theory (Figure 12).The condensate particle radii may be in the range 10 −2 − 10 −1 µm, and the pressure at the clouds may be in the range 0.08 − 0.2 bar (Figure 13).-We demonstrated that analytic phase curves with isotropic multiple scattering are in excellent agreement with full Rayleigh multiple scattering calculations (Figure 9).

Fig. 1 .
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 A g = 1 (upper limit).

Fig. 3 .
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 C 11 and A g , 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 x 1 , x 2 , and g.

Fig. 4 .
Fig.4.Kepler 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-minute 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.Article number, page 7 of 22

Fig. 5 .
Fig.5.Kepler SAP fluxes (gray points) after CBV corrections are applied, in relative flux units of ppm, for Quarter (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.2.Bayesian model selectionWhen 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.
3.1), we can efficiently compute the posterior predictive density for the full model including the GP, and determine Parameter TrES-2 b Kepler-7 b Kepler-41 b HAT-P-7 b KOI-13 b Model comparison results

Fig. 6 .Fig. 7 .
Fig.6.Inferred reflectance and temperature maps for each planet, where the center of the map corresponds to the sub-stellar 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 x 1 and x 2 .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. 8 .
Fig. 8. Stability curves of selected high-temperature condensates for an atmosphere with the metallicity of Kepler-41 b [M/H] = 0.38 (Bonomoet al. 2015).The calculations are assuming chemical equilibrium in the gas phase.

Fig. 10 .
Fig. 10.Geometric albedos for four Kepler targets, measured in four one-year segments of the Kepler light curves.

Fig. 11 .
Fig. 11.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 Figure12.
Fig.12.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 2018b, 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 Figure11) 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 CaTiO 3 (perovskite), Mg 2 SiO 4 (forsterite), and MgSiO 3 (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 Figure13.

Fig. 13 .Fig. 14 .
Fig. 13.Joint posterior distributions for the planetary longitudes bounding the less reflective region (x 1 , x 2 ) 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 (x 1 , x 2 ).Each species' condensation temperature is converted to a coordinate in (x 1 , x 2 ), 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 CaTiO 3 clouds located near 0.19 bar.Plausible cloud pressures for all species range from 0.08 to 0.2 bar.
found that the thermal flux of HD 189733 b and HD 209458 b change by less Article number, page 16 of 22 Brett M. Morris et al.: Observations of scattered light from exoplanet atmospheres

Table 1 .
Prior Morris et al. (2022)ues adopted for the reflected light fromHeng et al. (2021)and thermal emission phase curves fromMorris et al. (2022), as well as parameters for stellar and instrumental artifacts.⋆ The inhomogeneous reflector model defines A g as a free parameter, whereas the homogeneous reflector model derives a deterministic A g from ω, g, and the scattering phase function.
Betancourt 2017) 2018;Phan et al. 2019)nd 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.numpyroprovidessupportfor 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, seeBetancourt 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.4 Article number, page 4 of 22 Brett M.Morris et al.:Observations of scattered light from exoplanet atmospheres 2.9.Posterior Sampling Article number, page 12 of 22 Brett M. Morris et al.: Observations of scattered light from exoplanet atmospheres