# Spatially Resolved Modeling of Optical Albedos for a Sample of Six Hot
Jupiters

Danica J. Adams Division of Geological and Planetary Sciences, California
Institute of Technology, Pasadena, CA 91125, USA  Tiffany Kataria Jet
Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109,
USA  Natasha E. Batalha NASA Ames Research Center, MS 245-3, Moffett Field, CA
94035, USA  Peter Gao Department of Astronomy and Astrophysics, University of
California, Santa Cruz, Santa Cruz, CA 95064, USA  Earth and Planets
Laboratory, Carnegie Institution for Science, 5241 Broad Branch Road, NW,
Washington, DC 20015, USA  Heather A. Knutson Division of Geological and
Planetary Sciences, California Institute of Technology, Pasadena, CA 91125,
USA

###### Abstract

Optical secondary eclipse measurements made by _Kepler_ reveal a diverse set
of geometric albedos for hot Jupiters with equilibrium temperatures between
$1550-1700$ K. The presence or absence of high altitude condensates, such as
Mg2SiO4, Fe, Al2O3, and TiO2, can significantly alter optical albedos, but
these clouds are expected to be confined to localized regions in the
atmospheres of these tidally locked planets. Here, we present 3D general
circulation models and corresponding cloud and albedo maps for six hot
Jupiters with measured optical albedos in this temperature range. We find that
the observed optical albedos of K2-31b and K2-107b are best matched by either
cloud free models or models with relatively compact cloud layers, while
Kepler-8b and Kepler-17b’s optical albedos can be matched by moderately
extended ($f_{sed}$ = 0.1) parametric cloud models. HATS-11b has a high
optical albedo, corresponding to models with bright Mg2SiO4 clouds extending
to very low pressures ($f_{sed}$ = 0.03). We are unable to reproduce
Kepler-7b’s high albedo, as our models predict that the dayside will be
dominated by dark Al2O3 clouds at most longitudes. We compare our parametric
cloud model with a two-zone microphysical cloud model (CARMA). We find that
even after accounting for the 3D thermal structure, no single cloud model can
explain the full range of observed albedos within the sample. We conclude that
a better knowledge of the vertical mixing profiles, cloud radiative feedback,
cloud condensate properties, and atmospheric metallicities is needed in order
to explain the unexpected diversity of albedos in this temperature range.

keywords

††journal: ApJ††software: numba (Lam et al., 2015), pandas (McKinney, 2010),
bokeh (Bokeh Development Team, 2014), NumPy (Walt et al., 2011), IPython
(Pérez & Granger, 2007), Jupyter, (Kluyver et al., 2016), Virga (Batalha et
al., 2020), PICASO (Batalha et al., 2019)

##  1 Introduction

Transiting short period gas giant planets, or ‘hot Jupiters’, are among the
most favorable targets for atmospheric characterization studies. By observing
the decrease in optical flux when the planet passes behind the host star
(secondary eclipse), we can measure the optical dayside albedos for these
tidally locked planets. Most of these albedo measurements have come from space
telescopes observing in broad optical bandpasses (e.g. Rowe et al., 2008;
Kipping & Bakos, 2011; Demory et al., 2011, 2013; Coughlin & Lopez-Morales,
2012; Parviainen et al., 2013; Angerhausen et al., 2014; von Paris et al.,
2016; Niraula et al., 2018). The hot Jupiter HD 189733b is currently the only
planet with a spectroscopically resolved reflected light measurement (Evans et
al., 2013, geometric albedo of $0.40\pm 0.12$ at $290-450$ nm and $<0.12$ at
$450-570$ nm). WASP-12b and WASP-43b were also observed with optical
spectrographs on _HST_ , but both observations resulted in non-detections
(Bell et al., 2017; Fraine et al., 2021, geometric albedos $<0.06$). These
spectroscopic albedo measurements are invaluable for constraining the nature
of the scattering particles in the atmospheres of these planets (e.g. Barstow
et al., 2014).

Theoretical models predict that variations in hot Jupiter optical albedos
should primarily be driven by the presence or absence of high altitude
aerosols, which are expected to scatter incident starlight (e.g. Seager et
al., 2000; Burrows et al., 2008). Cloud-free atmosphere models for hot
Jupiters generally predict dayside albedos of less than 0.1 (Seager et al.,
2000), but the observed optical geometric albedo measurements published to
date span a wide range of values, with the brightest planets exhibiting
albedos as high as $0.3$ (e.g. Heng & Demory, 2013; Niraula et al., 2018). In
the solar system, the presence of water clouds increases Earth’s geometric
albedo to approximately 0.37 (e.g. Goode et al., 2001), while ammonia clouds
in Jupiter’s atmosphere contribute to its geometric albedo of approximately
0.5 (e.g. Marley et al., 1999). For the same reason, models suggest that the
presence of reflective condensates such as Mg2SiO4 or MgSiO3 in hot Jupiter
atmospheres can increase their albedos to values as high as 0.5 (Parmentier et
al., 2016, 2021; Roman et al., 2021). Unlike brown dwarfs, whose cloud
distributions and optical depths correlate closely with their equilibrium
temperatures (e.g. Kirkpatrick, 2005; Marley et al., 2010), hot Jupiter
albedos can vary by as much as an order of magnitude within a relatively
narrow range of equilibrium temperatures (Figure 1). This suggests that
individual planets with similar equilibrium temperatures may exhibit diverse
cloud properties.

Figure 1: Optical geometric albedo measurements for a sample of hot Jupiters
observed in the _Kepler_ bandpass; at these temperatures thermal emission is
negligible for planets with cloudy skies and the measured secondary eclipse
depth is dominated by reflected light. The equilibrium temperature is
calculated assuming a Bond albedo of zero and efficient day-night circulation.
The planet colors vary as a function of surface gravity. The dashed box
indicates the temperature range considered in this study. Measurements are
drawn from Fortney et al. (2011); Desert et al. (2011); Barclay et al. (2012);
Esteves et al. (2013); Shporer & Hu (2015); Niraula et al. (2018); Heng et al.
(2021).

Hot Jupiters are expected to be tidally locked as a result of their short
orbital periods, which can lead to significant day-night temperature gradients
(e.g. Guillot & Showman, 2002). This inhomogneous temperature structure
affects the cloud distribution, and 3D atmospheric circulation models predict
that hot Jupiters should host spatially inhomogenous clouds over a wide range
of equilibrium temperatures (e.g. Parmentier et al., 2016, 2021; Roman et al.,
2021). Furthermore, there is direct observational evidence for inhomogeneous
cloud structures stemming from reflected-light phase curves (e.g., Demory et
al., 2013). To date, three of the four planets with observed reflected-light
phase curves appear to have patchy clouds (Desert et al., 2011; Demory et al.,
2013; Angerhausen et al., 2014; Esteves et al., 2015; Shporer & Hu, 2015; von
Paris et al., 2016; Niraula et al., 2018, also see Figure 1). Kepler-7b, -12b,
and -41b all appear to have patchy clouds that are preferentially concentrated
in the western (dawn) hemisphere, causing the peak of the phase curve to occur
after the secondary eclipse (Shporer & Hu, 2015). Observations of the fourth
planet, TrES-2b, indicate that it is uniformly dark ($A_{s}<0.03$, which is
equivalent to $A_{g}<0.02$ for a Lambertian sphere, where $A_{s}$ is the
spherical albedo and $A_{g}$ is the geometric albedo) at all phases (von Paris
et al., 2016). These studies demonstrate that if we wish to explain the
observed planet-to-planet variations in the measured dayside optical albedos
of hot Jupiters, we must utilize 3D models capable of capturing the spatially
varying cloud structure.

While an increasing number of studies are using general circulation models
(GCMs) to predict cloud patterns, these have either focused on individual
planets (e.g. Webber et al., 2015; Lee et al., 2015; Oreshenko et al., 2016;
Helling et al., 2019; Helling, 2019; Helling et al., 2020; Lines et al., 2018,
2019; Roman & Rauscher, 2019) or generic grids of models (e.g. Parmentier et
al., 2016, 2021; Roman et al., 2021). Of these listed studies, only five
(Webber et al., 2015; Oreshenko et al., 2016; Parmentier et al., 2016, 2021;
Roman et al., 2021) use their models to calculate predicted albedos in the
_Kepler_ bandpass. The question of why individual planets with similar
equilibrium temperatures would exhibit widely varying cloud properties
therefore remains largely unexplored in the literature.

Table 1: Properties of our planet sample.

Planet |  $M_{P}$ ($M_{J}$) |  $R_{P}$ ($R_{J}$) |  ${T_{eq}}(K)$aaCalculated assuming a Bond albedo of zero and efficient day-night recirculation. |  $g$ $(m/s^{2})$ | a (AU) | Period (days) | Measured $A_{g}$ | Ref.1
---|---|---|---|---|---|---|---|---
K2-31b | 1.77 | 1.06 | 1550 | 39.07 | 0.022 | 1.26 | 0.023$\pm$0.002 | 1,9
Kepler-17b | 2.45 | 1.31 | 1570 | 35.41 | 0.026 | 1.49 | 0.099$\pm$0.017 | 2,3
HATS-11b | 0.85 | 1.51 | 1560 | 7.97 | 0.051 | 3.62 | 0.270$\pm$0.052 | 4,5,9
Kepler-7b | 0.45 | 1.65 | 1630 | 4.16 | 0.062 | 4.89 | 0.25$\pm$0.01 | 6,8
K2-107b | 0.84 | 1.44 | 1650 | 10.26 | 0.048 | 3.31 | 0.062$\pm$0.010 | 5,7,9
Kepler-8b | 0.60 | 1.41 | 1680 | 7.32 | 0.047 | 3.85 | 0.124$\pm$0.013 | 6,10

Note. — (1) Grziwa et al. (2016), (2) Desert et al. (2011), (3) Bonomo, A. S.
et al. (2017), (4) Bayliss et al. (2018), (5) Livingston et al. (2018), (6)
Esteves et al. (2015), (7) Eigmuller et al. (2017), (8) Heng et al. (2021),
(9) Niraula et al. (2018), (10) Esteves et al. (2013)

Modeling the 3D structure of clouds can be computationally demanding,
especially given the large number of parameters that must be considered
(including but not limited to the particle number density and size
distribution, spatial extent of the clouds, location of the cloud decks,
number and composition of cloud species, and coupled radiative feedback). In
examining the body of published studies that use GCMs to predict hot Jupiter
dayside albedos, we find that all of these studies make simplifying
assumptions in order to conserve run time and numerical complexity. Oreshenko
et al. (2016), for example, determined the locations of clouds by comparing
the 3D thermal structure from a GCM with relevant condensation curves, but
assumed a fixed particle size and local condensation only (i.e., they
neglected vertical mixing) for the cloud layers. In Parmentier et al. (2016),
cloud distributions were determined using the completed temperature structure
calculated by a non-gray cloud-free GCM, with a cloud-top minimum pressure
fixed at 1 microbar. In Parmentier et al. (2021), cloud layers were also
calculated using the thermal structure output from a GCM, but the prescribed
vertical extent is limited by temperatures within the extent of 200 mbar to 1
$\mu$bar. Roman & Rauscher (2019) and Roman et al. (2021) differ significantly
from the Parmentier models. In the former, cloud distributions were determined
at each timestep in a double-grey GCM which included radiative feedback. They
included extended clouds which form when the temperature-pressure profile
permits, but are forced to taper off at pressures between 0.3 mbar and 0.057
mbar. They also consider compact cases with varied optical thickness where the
clouds are truncated and tapered off after approximately one scale height
regardless of where their base forms. In Roman et al. (2021), they also
consider a more extensive grid of planet models, while also varying cloud
compositions, densities, and vertical extents. All of the studies listed above
assumed homogeneous condensation, allowing them to treat each individual cloud
species separately. Although there are studies in the literature that have
combined GCMs with microphysical cloud models (e.g., Lee et al., 2015; Lines
et al., 2018; Helling et al., 2019; Helling, 2019; Helling et al., 2020),
these models focused on individual planets and were limited in their albedo
predictions in the _Kepler_ bandpass. Although there are studies in the
literature that have combined GCMs with microphysical cloud models (e.g., Lee
et al., 2015; Lines et al., 2018; Helling et al., 2019; Helling, 2019; Helling
et al., 2020), these models focused on individual planets and were limited in
their albedo predictions in the _Kepler_ bandpass.

In this study, we utilize a suite of models to investigate the role of patchy
clouds over a sample of six individual hot Jupiters, chosen due to their
diverse observed albedos over a narrow range of equilibrium temperatures. In
§2, we describe our sample selection and summarize our modeling approach. We
use 3D GCMs to derive the thermal structure and eddy diffusion coefficients
($K_{zz}$). We then use Virga, a phase equilibrium cloud code, to make
detailed maps of cloud structure over the dayside of each planet. We use
PICASO, a radiative transfer program, to calculate the corresponding geometric
albedo maps and hemisphere-integrated dayside albedos. We also consider a
microphysical code, CARMA, which enables us to take a closer look at the role
of nucleation, condensation, and sedimentation in shaping the distribution of
dayside clouds. In §3 we compare our predicted dayside-integrated optical
albedos with observations of the six planets of interest. We then investigate
the relative importance of various model assumptions, such as equilibrium
condensation versus kinetic condensation, by comparing the Virga results with
the CARMA results. Finally, we discuss the implications of our results in §4
and present our conclusions in §5.

##  2 Methods

Figure 2: Map of the dayside temperatures (in Kelvin) of each planet at 1
mbar, roughly the pressure of unit optical depth in a clear atmosphere in the
_Kepler_ bandpass. Each planet is given a unique scale for the color bar to
best match the relevant temperature range.

###  2.1 Planet Sample

In this study we focus on planets with equilibrium temperatures ($T_{eq}$,
calculated assuming an albedo of zero and efficient day-night recirculation)
between $1550-1700$ K. For planets in this relatively narrow temperature
range, we expect that reflective silicate clouds should dominate the optical
dayside albedos (e.g. Parmentier et al., 2016; Powell et al., 2018; Parmentier
et al., 2021; Gao et al., 2020; Roman et al., 2021). This temperature range
contains some of the most reflective hot Jupiters observed to date, including
Kepler-7b (Demory et al., 2013; Heng et al., 2021) and HATS-11b (Niraula et
al., 2018). Note that two observed geometric albedos have been reported for
Kepler-7b, and from here on we consider the most recent value from Heng et al.
(2021). It also includes three moderately reflective hot Jupiters (Kepler-8b,
Kepler-17b, and K2-107b; Esteves et al., 2013; Desert et al., 2011; Niraula et
al., 2018) and one very dark hot Jupiter (K2-31b; Niraula et al., 2018).
Previous studies (e.g., Demory et al., 2013) have concluded that planets with
these equilibrium temperatures should have negligible amounts of thermal
emission in the _Kepler_ bandpass. We checked this by using PICASO to compute
the predicted thermal emission for each of the six planets in our sample and
found a contribution of 2 ppm or less. This is significantly smaller than the
uncertainties on the secondary eclipse depths used to calculate the geometric
albedos for these planets. Hence, throughout this paper we ignore the thermal
contribution to the measured secondary eclipse depth in the _Kepler_ bandpass.

In this study we present a GCM tailored to each of the six individual planets,
which we use to predict planet-to-planet variations in cloud coverage and
dayside albedo in the optical _Kepler_ band. These planets sample a range of
surface gravities, allowing us to investigate its effect on the planet’s
thermal structure, cloud distributions, and optical albedo. We summarize the
physical and orbital properties of each system in Table 1.

###  2.2 Modeling Atmospheric Circulation with the SPARC/MITgcm

We model each planet’s clear-sky (cloud-free) 3D thermal structure and
atmospheric circulation with the Substellar and Planetary Radiation and
Circulation (SPARC) model, which couples the GCM maintained at the
Massachusetts Institute of Technology (the MITgcm; Adcrot et al., 2004) with a
plane-parallel, two-stream version of the multi-stream radiation code as
described in Marley et al. (1999). The MITgcm solves the 3D primitive
equations on a staggered Arakawa C grid (Arakawa et al., 1977) with the
finite-volume method. The equations are discretized on a 128 $\times$ 64
cubed-sphere grid with 53 vertical layers extending from 200 bars at the
bottom boundary to 20 $\mu$bar at the top boundary. A horizontal fourth-order
Shapiro filter is used to smooth horizontal noise. We let each model run for a
simulated 1000+ Earth days so as to reach quasi-steady state equilibrium.

The radiative transfer scheme solves the two-stream radiative transfer
equations using the correlated-k method (Goody et al., 1989; Marley et al.,
1999) over 11 spectral bins (Kataria et al., 2013). This coupling allows for
the self-consistent calculation of the heating and cooling rates of the
atmosphere with latitude, longitude and pressure. At each grid point, the
radiative transfer scheme calculates the upward and downward fluxes at each
pressure layer, which are used to update the heating/cooling rates. These
rates are used by the MITgcm to update the wind and temperature fields.
Opacities are computed at each pressure-temperature point assuming chemical
and thermodynamic equilibrium, using the solar photospheric elemental
abundances of Lodders (2003). We interpolate across the PHOENIX stellar
atmosphere models to generate an input spectrum for each host star. The
SPARC/MITgcm has been successfully utilized for a series of hot Jupiter
studies (e.g. Showman et al., 2009; Showman & Kaspi, 2013; Showman et al.,
2015; Kataria et al., 2013, 2015, 2016; Parmentier et al., 2013; Lewis et al.,
2014), and we refer the reader to Kataria et al. (2016) for further details.
Figure 2 shows the resulting thermal structure for each of the six planets in
our sample at 1 mbar (approximately the level of unit optical depth for clear
skies).

###  2.3 Computing Equilibrium Condensate Clouds with Virga

We use the thermal structure and vertical mixing rates from the SPARC/MITgcm
model as inputs to Virga (Batalha et al., 2020)111Code and documentation
available at https://natashabatalha.github.io/virga/, an open-source code that
calculates phase equilibrium cloud distributions. The cloud parameterization
used in this code is described in Ackerman & Marley (2001) and has been used
for 1D model studies across a wide range of exoplanets and brown dwarfs (e.g.
Fortney et al., 2006; Marley et al., 2010; Morley et al., 2015). This
parametric approach also allows us to sample the 3D cloud structure at a much
higher spatial resolution than for our microphysical models.

We calculate $K_{zz}$ profiles from the root-mean-square (rms) vertical
velocities derived from global horizontal averages at a given pressure level
from the GCMs by assuming $K_{zz}=w(z)L(z)$, where $w(z)$ is the horizontally
averaged global rms vertical velocity from the GCM simulations, $L(z)$ is
approximated as the atmospheric pressure scale height $H(z)$ (but could be a
fraction of $H(z)$; see Smith, 1998), and $z$ is altitude. Moses et al. (2011)
note that this is only an estimate; a better approach would involve
calculating $K_{zz}$ from the eddy vertical velocity times the eddy
displacement, but this information is not readily obtainable from the GCMs;
this could be resolved by adding passive tracers to future GCM models. Our
treatment may overestimate $K_{zz}$ in the $\sim$10–200 bar radiative region,
where the vertical motion often consists of small-scale wave oscillations.

We reduce the spatial resolution of our longitude and latitude grid from
128$\times$65 to 10$\times$10 by binning the pressure-temperature profiles and
corresponding Kzz profiles prior to running Virga. We retain the original 53
pressure levels in the rebinned grid. This binning has a negligible effect on
our calculation of the phase-integrated albedo and significantly reduces
computation time. We bin using the area mean with angles from the Chebyshev-
Gauss integration method that vary as a function of planetary latitude and
longitude.

In Virga the molar mixing ratio of the condensed phase, $q_{c}$, is calculated
by solving the equation:

$-K_{zz}\frac{\partial q_{t}}{\partial z}-f_{sed}\omega_{*}q_{c}=0$ (1)

where $q_{t}$ is the total mixing ratio (condensed and vapor phases),
$\omega_{*}$ is the convective velocity scale, and $f_{sed}$ is defined as the
ratio of the mass-weighted droplet sedimentation velocity to the convective
velocity, $\omega_{*}$. The product $f_{sed}\omega_{*}$ describes an average
sedimentation velocity for the condensate, which offsets turbulent mixing. We
refer the reader to Ackerman & Marley (2001) for more details regarding the
equations that govern Virga.

In Eq. 1, $f_{sed}$ is the only parameter that cannot be calculated directly
from the models. We therefore treat it as a tunable parameter and explore a
range of possible values. Models with larger values of $f_{sed}$ will have
high rates of sedimentation, concentrating the condensing species in the lower
atmosphere. Conversely, models with smaller values of $f_{sed}$ will have much
slower sedimentation rates, allowing cloud particles to remain lofted higher
in the atmosphere. For each planet, we run a suite of Virga models with
$f_{sed}$ values of 0.03, 0.1, 0.3, 1.0, 3.0, and 6.0. This range is motivated
by comparisons to observational data, which suggest that $f_{sed}$ can be as
small as 0.01 for super Earths (Morley et al., 2015) or as large as 2-5 for
some gas giants and brown dwarfs (Skemer et al., 2016; MacDonald et al., 2018;
Saumon & Marley, 2008). For Jupiter’s ammonia clouds, an $f_{sed}$ of $\sim$2
appears to provide the best match to observations (Ackerman & Marley, 2001).
We therefore conclude that our chosen list of values spans a representative
range for this parameter.

###  2.4 Computing Microphysical Clouds with CARMA

Figure 3: Condensation curves (dashed grey) of Fe, Al2O3, TiO2, and Mg2SiO4
clouds compared with the two-zone model averaged temperature profiles (see
§3.2). Dayside (solid) and western limb (dashed) are shown for each planet:
K2-107b (yellow), Kepler-8b (orange), HATS-11b (black), K2-31b (navy),
Kepler-17b (purple), and Kepler-7b (green).

In addition to the parametric model described above, we also utilize the
Community Aerosol and Radiation Model for Atmospheres (CARMA), a more
computationally demanding microphysical cloud formation model. CARMA
calculates the equilibrium cloud particle size distribution by solving the 1D
discretized continuity equation for aerosol particles that experience vertical
transport due to sedimentation and eddy diffusion and production and loss due
to particle nucleation (homogeneous and heterogeneous), condensation,
evaporation, and coagulation. CARMA has been previously used to investigate
condensate cloud formation on Earth (e.g., Ackerman et al., 1993, 1995; Jensen
& Toon, 1994), Venus (e.g. James et al., 1997; McGouldrick & Toon, 2007; Gao
et al., 2014), Mars (e.g. Colaprete et al., 1999), and exoplanets (e.g., Gao
et al., 2018, 2020; Powell et al., 2019). CARMA has also been used to model
photochemical hazes on Titan (Toon et al., 1992), Pluto (Gao et al., 2017),
ancient Earth (Wolf & Toon, 2010), and warm Jupiters (Adams et al., 2019). In
order to make our use of CARMA computationally tractable, we divide each
planet into two zones and calculate averaged temperature and $K_{zz}$ profiles
for each zone as described in §3.2 (see also Figure 3). We do not consider
photochemical hazes here, as the planets in our sample lie above the
temperature range where these hazes are expected to form (Lodders & Fegley,
2002; Kawashima & Ikoma, 2019; Gao et al., 2020).

In the CARMA model, the formation of condensate clouds begins with either
homogeneous or heterogeneous nucleation. Cloud particles undergo homogeneous
nucleation when stable clusters of condensate molecules form and grow directly
from the vapor. The rate is controlled by the material properties of the
condensate, such as its molecular weight and surface energy, and the flux of
molecules to the cluster, which depends on the abundance of condensate vapor.
Under the same supersaturation and local temperature, high surface energy and
molecular weight materials tend to nucleate more slowly than low surface
energy and low molecular weight materials. Unlike homogeneous nucleation,
heterogeneous nucleation involves the formation of stable clusters on
condensation nuclei, or foreign surfaces, which are provided by other aerosol
particles in the atmosphere. The size and abundance of these particles
strongly impact the rate of heterogeneous nucleation. The nucleation rate is
also dependent on the interaction between the condensate and the surface,
characterized by the contact angle between the surface and the condensate
cluster, the energy needed by a condensate molecule to desorb from the
surface, and the oscillation frequency of the condensate molecule on the
surface, which is related to the desorption energy (Pruppacher & Klett, 1978).

Unlike in Virga, which assumes a log-normal particle size distribution,
particle size distributions in CARMA are resolved using mass bins and can
change over time via condensation, evaporation, and coagulation. The former
two rates are dependent on the flux of condensate molecules and the rate at
which particles may conduct away the latent heat released upon condensation.
Coagulation, or growth via physically sticking upon the collision of
particles, is also considered. Each mass bin corresponds to particle masses
twice that of the previous bin. We use 65 bins in our model, with the mass in
the first bin corresponding to particles with radii of 0.1 nm for all species.

###  2.5 Cloud Compositions

Figure 4: Nadir optical depths (integrated over the Kepler bandpass) for three
condensate species in our Virga model of K2-107b. From left to right: Mg2SiO4,
TiO2, Al2O3. Each panel shows two representative grid points: black
corresponds to a western grid point centered at $-60\arcdeg$ W, $8\arcdeg$ N
and red corresponds to a dayside gridpoint centered at $42\arcdeg$ E,
$8\arcdeg$ N. The value of $f_{sed}$, which ranges from 0.03 to 6.0, is
indicated by the shading of each line, where the lightest shades correspond to
the smallest values and darkest shades represent the largest values of
$f_{sed}$.

For the Virga cloud modeling in this study, we are only interested in cloud
species that are expected to be abundant in these atmospheres and which
condense at relatively low pressures (approximately 1 bar). We identify three
cloud species that are likely to be important based on comparing the species’
condensation curves to the planets’ temperature pressure profiles (Figure 3;
also see Ackerman & Marley, 2001; Morley et al., 2012): Mg2SiO4, Al2O3, and
TiO2. This assumes that all SiO goes into Mg2SiO4 rather than SiO2 or MgSiO3.
The inclusion of SiO2 would require a different modeling framework including
kinetic condensation. We plot the optical depths in the Kepler bandpass (shown
in Figure 6) for these species as a function of pressure as computed by Virga
at a representative location on the planet K2-107b in Figure 4. Although some
previous studies of spatially inhomogeneous cloud formation in hot Jupiters
also included sulfide clouds (e.g. Parmentier et al., 2016), recent
microphysical cloud models (Powell et al., 2019; Gao et al., 2020) indicate
that the high surface energy of sulfide condensates makes homogeneous
nucleation unlikely. This conclusion is in good agreement with observational
data indicating that the observed cloud opacity for more moderately irradiated
hot Jupiters is well-matched by models without sulfide condensates (e.g.
Chachan et al., 2019; Gao et al., 2020). We therefore exclude sulfide
condensates from our Virga models. Fe is also a potentially important
condensate, but like the sulfides previous microphysical modeling suggested
that Fe cloud formation proceeds slowly due to Fe’s high surface energy (Gao
et al., 2020). On the other hand, Fe has been considered a major cloud
component in many previous works (Marley et al., 2013; Helling et al., 2008;
Marley & Robinson, 2015; Fortney, 2005; Sudarsky et al., 2003). Therefore,
while our nominal Virga models will not include Fe, we will consider another
set of simulations where Fe is included to explore its impact on our predicted
albedos.

For the CARMA simulations in this study, we use the same model setup as in Gao
et al. (2020) where clouds of TiO2, Al2O3, Fe, Mg2SiO4, Cr, MnS, Na2S, ZnS,
and KCl are allowed to form. However, as previously discussed, in practice
only Mg2SiO4, Al2O3, and TiO2 form in any significant abundance. Of these
three species, Mg2SiO4 is particularly unlikely to undergo homogeneous
nucleation, as it is not abundant in the vapor phase. Instead, it is the
product of a thermochemical reaction between Mg, SiO, and H2O (e.g. Visscher
et al., 2010). Similarly, Al2O3 does not exist in the gaseous phase as other
aluminum oxide species will be more stable (e.g. Patzer et al., 2005). As in
Gao et al. (2020), we allow these two species to heterogeneously nucleate on
homogeneously nucleated TiO2 seeds. Though Fe is permitted to both
homogeneously and heterogeneously nucleate (on TiO2 seeds), the high surface
of Fe prevents significant Fe cloud formation.

###  2.6 Computing Optical Albedos with PICASO

(a)

(b)

Figure 5: Pressure level (bar) of unit optical depth for Rayleigh scattering
(pink), cloud opacity (green), and molecular opacity (black) as a function of
wavelength for a single grid point ($-60\arcdeg$ W, $8\arcdeg$ N) in the Virga
model of hot Jupiter K2-107b. $f_{sed}$ of 0.1 is shown on the left and 3.0 on
the right. The shaded regions indicate the dominant opacity source as a
function of wavelength. This is the same grid point shown in Figure 4.

We use the outputs from the cloud codes Virga and CARMA to calculate the
single scattering albedos, asymmetry parameters, and optical depths for each
pressure layer at each location in the atmosphere assuming Mie scattering. We
then convert these quantities into wavelength-dependent albedos using The
Planetary Intensity Code for Atmospheric Scattering Observations (PICASO;
Batalha et al., 2019). This code is governed by the radiative transfer
equation:

$\begin{gathered}I(\tau_{i},\mu)=I(\tau_{i+1},\mu)\exp(\frac{\delta\tau_{i}}{\mu})-\\\
\int_{0}^{\delta\tau_{i}}S(\tau^{\prime}\mu)\exp(-\frac{\tau}{\mu})d\tau^{\prime}/\mu\end{gathered}$
(2)

where $I(\tau_{i},\mu)$ is the azimuthally averaged intensity emerging from
the top of a layer $i$ with opacity $\tau_{i}$ and outgoing angle $\mu$.
$I(\tau_{i+1})\exp(\frac{\delta\tau_{i}}{\mu})$ is the incident intensity at
the lower boundary of the layer attenuated by the optical depth within the
layer $\delta\tau$, and $S(\tau^{\prime},\mu)$ is the source function
integrated over all layers. The source function has two components: single
scattered and multiple scattered radiation integrated over all diffuse angles:

$\begin{gathered}S(\tau^{\prime},\mu)=\frac{\omega}{4\pi}F_{o}P_{single}(\mu,-\mu_{o})\exp(-\frac{\tau^{\prime}}{\mu_{s}})+\\\
\frac{\omega}{2}\int_{-1}^{1}I(\tau^{\prime},\mu^{\prime})P_{multi}(\mu,\mu^{\prime})d\mu^{\prime},\end{gathered}$
(3)

where $\omega$ is the single scattering albedo and $F_{o}$ is the incident
flux. $P_{multi}$ and $P_{single}$ describe the phase function of the multiple
and single scattering, respectively. $P_{single}$ is an opacity weighted
combination of the Rayleigh phase function and a two-term Henyey Greenstein
phase function. $P_{multi}$ requires integration over all diffuse angles, for
which PICASO uses a $N=2$ Legendre expansion. This approximation alone is
inadequate to represent cases with high rates of forward scattering, so PICASO
implements the delta-Eddington approximation to scale $g$, $\omega$, and
$\tau$ to more accurately capture the forward scattering peak.

PICASO considers the extinction from three opacity sources in order to
calculate the geometric albedo as a function of wavelength: molecular
absorption, Rayleigh scattering, and scattering by clouds. We show a
representative calculation of these three opacity sources as a function of
wavelength in Figure 5. To describe the phase-dependence, PICASO computes the
emergent intensity from the disk at multiple plane-parallel facets, where each
has its own incident and outgoing angles. PICASO uses the Chebyshev-Gauss
integration method to integrate over all emergent intensities. We integrate
the wavelength-dependent geometric albedo over the _Kepler_ response function
for each point in our $10\times 10$ grid and then integrate again over the
dayside hemisphere in order to obtain a geometric albedo that we can compare
with the _Kepler_ measurements. We have run sensitivity tests that demonstrate
a higher resolution grid ($20\times 20$ grid) yields comparable results.

###  2.7 Effect of Simplifying Model Assumptions

In this study we do not consider radiative feedback from clouds, which might
affect our albedo predictions. Clouds can alter the planet’s global thermal
structure in several ways. First, they can reduce the amount of heating on the
dayside by increasing the planet’s Bond albedo. Second, they can suppress
cooling on the nightside by preventing the re-radiation of infrared light to
space. Parmentier et al. (2021) and Roman et al. (2021) ran grids of 3D GCM
models incorporating radiative feedback from clouds spanning a range of
incident fluxes. They found that the presence of reflective clouds on the
dayside resulted in lower global temperatures, while the presence of nightside
clouds inhibited cooling, causing a $100-200$ K global increase in
temperature. With the possible exception of Kepler-7b (see §4.2.1), we expect
that such shifts in temperature would not substantially alter the pressures of
the cloud decks or reduce their horizontal extent for the planets examined
here.

We also note that in this study, our 3D atmospheric circulation models are
decoupled from our cloud models. While this does enable more flexibility in
exploring different cloud species and sedimentation parameters in both Virga
and CARMA, cloud formation and transport is ultimately a coupled process
between advection, radiation, and chemistry. Previous studies that couple
cloud microphysical models and atmospheric circulation models (e.g., Lee et
al., 2015; Lines et al., 2018) suggest that zonal (east/west) and/or
meridional (north-south) transport of cloud particles from colder regions of
the atmosphere could lead to high cloud particle number densities, even on the
hotter dayside, which could in turn enable more cloud nucleation and growth.
These processes will only affect our dayside albedo predictions if the planet
in question has large dayside temperature (and hence albedo) gradients with
longitude or latitude. We revisit both of these assumptions (radiative
feedback and cloud microphysics coupled to circulation) in §4.2, where we
discuss their implications for individual planets in light of our model
results.

##  3 Results

The SPARC/MITgcm models indicate that the day-night temperature contrasts for
the six planets in our sample vary in magnitude (Figure 2). As expected, the
hottest region on the day side in all six models is located to the east of the
substellar point. This shift is caused by eastward equatorial winds, which
transport heat to the planet’s nightside (e.g., Showman et al., 2020, and
references therein).

Kepler-7b has the largest thermal gradient of all the planets in our sample,
followed by K2-31b and Kepler-17b; this is because the former has a relatively
low surface gravity (approximately 4 m/s2) and the latter two have the
shortest orbital periods in the sample

Kepler-7b has the largest thermal gradient of all the planets in our sample,
followed by K2-31b and Kepler-17b; this is because the former has a relatively
low surface gravity (approximately 4 m/s2) and the latter two have the
shortest orbital periods in the sample. Our models indicate that Kepler-7b is
also warmer at depth than the other planets in this sample. This is expected,
as Kepler-7b has the lowest surface gravity of the six planets and under
hydrostatic equilibrium the gradient of temperature with respect to pressure
is related to the inverse of the surface gravity (e.g., Gao et al., 2018).

Figure 6: Hemisphere-averaged albedo of K2-107b as a function of wavelength
for a clear atmosphere (black) and cloudy atmospheres of varying $f_{sed}$
values (green; 0.03 as a solid line, 0.1 as a dashed line, and 0.3 as a dotted
line). Cloudy atmospheres including Fe condensates are shown in purple. All
cloud distributions are computed using Virga. The _Kepler_ response function
is overplotted in grey. Figure 7: Relative dayside albedo contribution at each
grid point including appropriate geometric weights. Planets are sorted by
increasing equilibrium temperature from top to bottom, and sorted by
increasing $f_{sed}$ (0.03, 0.1, 0.3, and 1.0) from left to right. We omit
$f_{sed}$ values of 3.0 and 6.0 from the figure, as the clouds in these models
reside below the level of unit molecular albedo, and thus the albedo remains
roughly unchanged compared to the $f_{sed}$ = 1.0 case. The relative albedo
contributions are all smaller than the face-integrated albedo by roughly a
factor of $~{}\pi$; summing over the points yields the integrated albedo. The
models include three cloud species: Mg2SiO4, Al2O3, and TiO2.

###  3.1 Virga Model Results

In order to determine the effect of the clouds on the albedo, we must first
calculate their vertical extent at each location in our model grid. In our
Virga equilibrium cloud models, the vertical extend of the clouds is
controlled by our choice of $f_{sed}$. Taking K2-107b as a representative
example, Figure 4 indicates that Mg2SiO4 clouds reach a unit optical depth
near 1 mbar for the lowest $f_{sed}$ value, 0.03, and near 50 mbar for
$f_{sed}=$ 0.1. These clouds will only contribute to the albedo at wavelengths
where they reach optical depth unity at lower pressures than Rayleigh
scattering or molecular opacity (Figure 5). This means that clouds will have a
greater effect at wavelengths where the combined molecular and Rayleigh
scattering opacity is lower. We show the effect of varying $f_{sed}$ on the
wavelength-dependent albedo of K2-107b in Figure 6. As $f_{sed}$ decreases and
the clouds extend to lower pressures, the cloud opacity contributes most to
the overall albedo for an increasingly large fraction of the _Kepler_
bandpass.

We next examine how the contribution of clouds to the band-integrated albedo
varies across the dayside atmosphere. In Figure 7, we show the $10\times 10$
grid of albedos in the _Kepler_ bandpass for each planet as a function of
$f_{sed}$. We find that the three planets with the greatest temperature
variation as a function of longitude (Kepler-7b, K2-31b and Kepler-17b) also
have relatively large albedo variations across their dayside atmospheres for
low to intermediate $f_{sed}$ values. Once $f_{sed}$ increases above 0.3, the
clouds remain confined below the optical depth unity level of molecular
absorbers, such that the presence or absence of clouds does not affect the
observed albedo. These three planets appear brighter on their western limbs
than in the east, in good agreement with the albedo maps derived from the
_Kepler_ phase curve for Kepler-7b (Demory et al., 2013) and other planets
with comparable equilibrium temperatures (Shporer & Hu, 2015). In contrast to
these three planets, HATS-11b, K2-107b, and Kepler-8b all appear to have
fairly homogeneous albedos, as expected based on their more homogeneous
thermal structures.

Table 2: Hemisphere-Averaged Albedos from Virga Models as a Function of $f_{sed}$aaBolded values indicate the simulated Virga albedo that best matches the _Kepler_ eclipse observations. Planet | Measured | Clear | 0.03 | 0.1 | 0.3 | 1.0 | 3.0 | 6.0 | 0.1, with Fe
---|---|---|---|---|---|---|---|---|---
K2-31b |  0.023$\pm$0.002 | 0.015 | 0.404 | 0.123 | 0.023 | 0.015 | 0.015 | 0.015 | 0.115
Kepler-17b |  0.099$\pm$0.017 | 0.017 | 0.416 | 0.131 | 0.023 | 0.017 | 0.017 | 0.017 | 0.125
HATS-11b |  0.270$\pm$0.052 | 0.066 | 0.301 | 0.127 | 0.060 | 0.066 | 0.066 | 0.066 | 0.127
Kepler-7b |  0.194$\pm$0.013 | 0.009 | 0.064 | 0.034 | 0.015 | 0.010 | 0.009 | 0.009 | 0.030
K2-107b |  0.062$\pm$0.010 | 0.065 | 0.333 | 0.154 | 0.063 | 0.065 | 0.065 | 0.065 | 0.134
Kepler-8b |  0.124$\pm$0.013 | 0.069 | 0.319 | 0.151 | 0.063 | 0.070 | 0.070 | 0.070 | 0.151
Figure 8: Geometric albedo predictions in the _Kepler_ optical bandpass
compared with published albedo measurements as a function of equilibrium
temperature. Albedo predictions from the full-resolution Virga models are
shown as circles where the shading indicates the $f_{sed}$ value, going from
0.03 (light) to 6.0 (dark). _Kepler_ albedo measurements are shown as stars
while the predicted albedos from the two-zone CARMA models are shown as
triangles. Left panel excludes Fe in Virga calculations, while the right panel
includes Fe condensates.

Lastly, we compare the hemisphere-averaged dayside albedo in the _Kepler_
bandpass as a function of $f_{sed}$ to the measured dayside albedo for each
planet (Figure 8). We find that K2-31b and K2-107b are best described by
models with large values of $f_{sed}$ or (equivalently) clear atmospheres,
indicating that any reflective clouds present in these atmospheres do not
extend above the level of unit molecular opacity. Kepler-17b and Kepler-8b are
best-matched by models with intermediate $f_{sed}$ values, implying that their
albedos are moderately enhanced by cloud opacity. HATS-11b is best matched by
models with low $f_{sed}$ values, corresponding to a high, vertically extended
reflective cloud layer spanning a wide range of longitudes. No $f_{sed}$ value
is able to reproduce the high observed albedo for Kepler-7b, which our models
indicate is too warm for bright Mg2SiO4 clouds to form over much of the day
side. Instead, our models indicate that the dayside of Kepler-7b is dominated
by deep Al2O3 clouds at pressures of around a bar. The addition of iron clouds
do not greatly change our results except for decreasing the albedos of low
$f_{sed}$ models, as iron clouds sink below the photosphere for higher
$f_{sed}$’s. We list the observed albedos and the predicted albedos as a
function of $f_{sed}$ in Table 2. We find that the hemisphere-integrated
albedo over the _Kepler_ bandpass is very sensitive to the assumed value of
$f_{sed}$.

###  3.2 Comparison to CARMA Microphysical Models

We find that dayside albedos can vary significantly depending on the assumed
sedimentation efficiency (e.g., Figure 7). We therefore use these maps to
divide each planet into two zones, and then utilize the more computationally
demanding CARMA microphysical cloud model to solve for the vertical extent of
the clouds and corresponding albedo in each zone.

Our Virga albedo maps indicate that a subset of the planets in our sample are
likely to have spatially inhomogenous Mg2SiO4 clouds located near their
western limbs. Although it would be computationally prohibitive to run a
separate CARMA model for each point in the $10\times 10$ grid, we can
nonetheless capture this cloud structure using a more computationally
tractable two-zone model. We define a threshold value in longitude for each
planet based on the albedo distributions found by Virga as shown in Figure 7
and in Table 3. We determine the longitude range defining the two zones (a
western zone and a dayside zone) as the division that yields the greatest
difference in albedo at mid-latitudes as determined from the Virga albedo
maps. Occasionally the western zone includes a limb region in the east; see
Table 3 for a list of the zone definitions for each planet. The resulting
zonally-averaged pressure-temperature profiles are plotted in Figure 3. We run
a separate CARMA model for each zone, and also run Virga models for the same
zones in order to facilitate comparisons between the two models. In the
discussion below, we limit our comparisons to these two-zone Virga models
unless otherwise noted.

Figure 9: Optical depth (left), single scattering albedo (center), and
asymmetry parameter (right) as a function of pressure for the two-zone Virga
(dashed; $f_{sed}$ fixed to 0.1) and CARMA (solid) models, calculated by
dividing the dayside hemisphere into a dayside (red) and western (black) zone
(see §3.2). Each row corresponds to a different planet, sorted left-to-right
by increasing equilibrium temperature. Figure 10: Dayside and western zone
albedos from the CARMA models integrated over the _Kepler_ bandpass. Virga
model albedos with $f_{sed}$ equal to 0.1 are listed below for comparison.
Both sets of models are calculated in two zonally averaged regions, defined in
Table 3. Planets are sorted from left to right and top to bottom by increasing
equilibrium temperature.

When comparing the two models, it is important to note that they predict
distinct particle compositions. While Virga assumes that particles nucleate
homogeneously, CARMA predicts that the Mg2SiO4 condensates will nucleate
heterogeneously onto TiO2 cores (see §2.5). We find that accounting for this
TiO2 core when we calculate the single scattering albedos for the Mg2SiO4
particles in our CARMA models results in a lower single scattering albedo than
when we repeat the calculation for pure Mg2SiO4 particles, which return single
scattering albedo profiles of greater than 0.96. For comparison, we refer the
reader to Figure 9; a single value for single scattering albedo is not
representative of heterogeneous particles, which are sensitive to the core
mass fraction of TiO2. However, even after accounting for this effect we find
that the species-averaged single scattering albedo shown in Figure 9 for our
CARMA models with heterogeneous particles is still comparable to the Virga
models with $f_{sed}$ equal to 0.1 (Figure 10).

When we calculate the hemisphere-integrated albedos for these planets from the
CARMA models we find that they generally lie within the range of Virga
predictions for varying $f_{sed}$ (see Figure 8 and Table 3). In Figure 8, we
show the hemisphere-integrated albedos from the full-resolution Virga model
from §2.3 rather than the two-zone model discussed in this section.

We find that K2-107b and HATS-11b appear to have the brightest hemisphere-
integrated CARMA albedos, with values consistent with those of the lowest
$f_{sed}$ ($0.03-0.1$) Virga models. Kepler-8b, Kepler-17b, and K2-31b have
somewhat lower predicted CARMA albedos, more comparable to the $f_{sed}$ of
$0.3$ Virga models. Kepler-7b is a notable exception, as our CARMA model
predicts an albedo that is higher than that of the smallest $f_{sed}$ Virga
model.

Table 3: Hemisphere-Averaged Dayside Albedos from CARMA Models

Planet | Measured | Clear | West | Day |  Day Zone11Reference to geometric albedo measurement is bolded for each row. | Integrated
---|---|---|---|---|---|---
K2-31b |  0.023$\pm$0.002 | 0.015 | 0.085 | 0.054 | $[-17\arcdeg W,51\arcdeg E]$ | 0.073
Kepler-17b |  0.099$\pm$0.017 | 0.017 | 0.074 | 0.069 | $[0\arcdeg E,90\arcdeg E]$ | 0.073
HATS-11B |  0.270$\pm$0.052 | 0.066 | 0.169 | 0.165 | $[51\arcdeg E,90\arcdeg E]$ | 0.168
Kepler-7b |  0.194$\pm$0.013 | 0.009 | 0.199 | 0.034 | $[-38\arcdeg W,51\arcdeg E]$ | 0.128
K2-107b |  0.062$\pm$0.010 | 0.065 | 0.233 | 0.229 | $[51\arcdeg E,90\arcdeg E]$ | 0.280
Kepler-8b |  0.124$\pm$0.013 | 0.069 | 0.082 | 0.109 | $[51\arcdeg E,90\arcdeg E]$ | 0.092

11footnotetext: Eastern/day zone boundaries are listed in the table. The
western/limb zone is defined as the region excluded by the dayside definition.

##  4 Discussion

###  4.1 Comparison to Previous Modeling Efforts

We can compare our results to generic grids of cloudy GCM models in the
published literature. First, we note that our GCM outputs for Kepler-7b
roughly agree with those shown in Oreshenko et al. (2016), and the general
distribution of our silicate and corundum clouds agrees with the modeling
results of Roman & Rauscher (2019).

Roman et al. (2021) investigated planetary albedos using a grid of GCMs with
varying irradiation temperature and surface gravity. The closest equivalent
models in their grid are for a planet with an irradiation temperature of
either 2500 K (equilibrium temperature of 1500 K) or 2750 K (equilibrium
temperature of 1700 K) and a surface gravity of 10 m/s2. We compare to their
nucleation-limited models, which also exclude iron and sulfide condensates.

The models presented in Roman et al. (2021) assume a fixed pressure-dependent
particle size for the clouds, with a size of 0.1 $\mu$m at the top of the
atmosphere that increases exponentially with increasing pressure for pressures
greater than 10 mbar. Since their equilibrium cloud models do not solve for
the vertical extent of the cloud layers, they present two cases corresponding
to compact (cloud tops limited to 1.4 scale heights above the cloud base) and
vertically extended (cloud top pressure of 0.1 mbar) cloud layers. Unlike our
models, they account for radiative feedback from these clouds when solving for
the temperature structure of the atmosphere. These models indicate that
compact cloud layers will result in relatively low and uniform dayside
albedos, in good agreement with our results (Figure 7). Their vertically
extended cloud models exhibit a range of optical spherical albedos between
$0.2-0.3$, with lower values for the higher temperature model as the
reflective silicate clouds become increasingly confined to the cooler western
region of the atmosphere. This also agrees with the qualitative picture from
our models, which span an equivalent range of albedos. They conclude that
their optical albedos for the nucleation-limited case are dominated by
silicate cloud particles, consistent with our conclusions here.

Parmentier et al. (2016) and Parmentier et al. (2021) considered GCMs spanning
a wide range of equilibrium temperatures. In Parmentier et al. (2016), the
clouds are post-processed (i.e., they do not include cloud radiative feedback
in the GCM models), and they model the condensation of a wide range of cloud
species. These models predict that MgSiO3 and CaTiO3 clouds should dominate
the dayside optical albedo for planets with equilibrium temperatures between
$1500-1700$ K. In Parmentier et al. (2021) they incorporate cloud radiative
feedback for the condensation of a single cloud species, MnS, but the
treatment of silicate clouds is similar to Parmentier et al. (2016). For that
reason, we focus on Parmentier et al. (2016) for our comparison. In this study
they assume a small fixed particle size distribution centered at 0.1 $\mu$m
and a cloud top pressure of 1 microbar, which yields geometric albedos greater
than 0.5 in the _Kepler_ bandpass for planets of approximately 1500 K
equilibrium temperature. This value is much higher than both our albedos and
those reported by Roman et al. (2021), and is most likely due to the very low
cloud top pressure assumed in these models. They invoke a cold trap for
silicates to reduce the albedo, while we predict that changing the cloud
vertical extent can achieve a similar difference in observable albedo.

###  4.2 Comparison to Published Kepler Albedos

Our results show that, with the exception of Kepler-7b, it is possible to
match the observed optical geometric albedos for all of the planets in our
sample using either Virga or CARMA models. However, no single model (Virga at
a fixed $f_{sed}$ or CARMA) can explain the observed albedos of all six
planets. Although the albedo predictions from the CARMA models are broadly
consistent with those of the Virga models, neither model is able to reproduce
or explain the observed planet-to-planet variations in dayside albedo. This
suggests that accounting for differences in equilibrium temperature, host star
spectral type, surface gravity, and rotation rate alone are not sufficient to
capture the observed diversity of hot Jupiter albedos in this temperature
range. We discuss our results in more detail on a planet-by-planet basis
below.

####  4.2.1 Kepler-7b: Spatially Resolved Clouds

Our CARMA and Virga models both indicate that Kepler-7b’s relatively bright
dayside albedo is dominated by reflection from the region near the western
limb, in good agreement with spatially resolved albedo constraints from phase
curve observations (Demory et al., 2013; Heng et al., 2021). This underscores
the importance of using spatially resolved cloud models for tidally locked hot
Jupiters. We note that there is some tension between our model predictions and
the observational data, as fits to Kepler-7b’s optical phase curve indicate
that the bright reflective western zone extends as far as $10\pm 6\arcdeg$
west of the substellar point (Muñoz & Isaak, 2015; Heng et al., 2021). Our
Virga models predict that the atmosphere will only be cool enough for Mg2SiO4
clouds to condense in the two westernmost longitude bins (extending from the
terminator to approximately 38$\arcdeg$ west of the substellar point; see
Figure 7 and Table 3). This is likely why our models under-predict Kepler-7b’s
optical geometric albedo.

We consider two possible explanations for this discrepancy. Zonal transport of
cloud particles from the western limb region could increase the albedo in
adjacent longitudes where the atmosphere is otherwise too warm for them to
condense (see §2.7). However, our models for Kepler-7b prefer small cloud
particles with a large vertical extent; these small particles might have a
relatively short lifetime in the hotter substellar region of the dayside
atmosphere. Instead, perhaps small particles transported meridionally could
nucleate and grow in bands at high latitudes (Lines et al., 2018). Ultimately,
this planet would be an interesting test case for microphysical transport
models (e.g., Lee et al., 2015; Lines et al., 2018), which can explicitly
quantify the timescales of these processes and predict the resulting
horizontal distribution of cloud particles.

Alternatively, if Kepler-7b’s dayside is cooler than predicted by our GCM,
Mg2SiO4 would be able to condense over a wider range of longitudes. Our GCMs
do not account for reflectivity from clouds when calculating the effect of
incident starlight on the dayside atmosphere; this effect might reduce the
magnitude of dayside heating and result in globally lower temperatures (Lines
et al., 2018; Roman & Rauscher, 2019; Roman et al., 2021). However, if the
clouds extend over a significant fraction of the planet’s nightside it could
result in net global warming, as they would act to reduce the amount of energy
that can be radiated to space in this region (Roman & Rauscher, 2019; Roman et
al., 2021; Parmentier et al., 2021). For Kepler-7b, whose clouds extend over
much of the western hemisphere, it is unclear which of these two competing
effects would dominate. These explanations assume that the dayside coverage
area is the most significant limiting factor on the brightness of the dayside-
integrated albedo. However, it is also important to consider factors that
might increase the brightness of the cloudy region, including a larger
vertical extent for the clouds or brighter cloud particles. We discuss how
porous particles may increase HATS-11b’s dayside albedo in Section 4.2.3; this
same explanation might also apply to Kepler-7b..

####  4.2.2 K2-31b and K2-107b Do Not Host High Altitude Reflective Cloud
Layers

The observed albedos of K2-31b and K2-107b are relatively low, and are well-
matched by Virga models with clear skies and/or deep clouds (i.e., those with
opacities dominated by molecular absorption). For K2-31b, the contribution of
clouds to the albedo is negligible for $f_{sed}$ of 0.3 and larger. For
K2-107b, we obtain a comparable result for $f_{sed}$ of 1.0 and larger. This
suggests that both of these planets may have relatively efficient
sedimentation (deep clouds), or alternatively that they have relatively little
condensable material in their atmospheres (perhaps corresponding to a
relatively low atmospheric metallicity). If the lack of clouds is due to
efficient sedimentation, this would appear to contradict predictions from our
microphysical CARMA models, which track the sedimentation of cloud particles
explicitly and predict albedos that are a factor of two or more higher than
the observed values for these two planets.

Our CARMA models utilize vertical mixing rates calculated from our GCMs. If
these mixing rates are overestimates of the true values, we might expect any
clouds near the day-night terminator on these planets to also be relatively
compact. If this is the case, the transmission spectra of these two planets
should show relatively strong absorption features. Although K2-31b has a high
surface gravity and is therefore a more challenging target for transmission
spectroscopy, K2-107b might be accessible to future space telescopes like
_JWST_. More broadly, the sedimentation rates calculated from GCMs and CARMA
could be tested with comparisons to transmission spectra from ongoing surveys
(Sing et al., 2016; Crossfield & Kreidberg, 2017; Fu et al., 2017). If CARMA
models underestimate the sedimentation efficiency for other planets, the
disagreement should be detectable in these data, which are very sensitive to
the vertical distribution of cloud particles near the limb. To date, most
planets appear to be well-matched by CARMA model predictions (e.g. Chachan et
al., 2020; Gao et al., 2020), indicating that such model-data disagreements
may be relatively rare.

####  4.2.3 HATS-11b, K2-31b, and K2-107b Have Reflective Dayside Clouds

The measured _Kepler_ albedo for HATS-11b ($0.27\pm 0.05$) is brighter than
that of Kepler-7b, in good agreement with our model predictions. This planet
is cooler than Kepler-7b, with a smaller day-night temperature gradient. As a
result, our models predict a global reflective dayside cloud layer for small
$f_{sed}$ values, and our hemisphere-integrated Virga albedo for the $f_{sed}$
of 0.03 model is within $1\sigma$ of the measured value. This is an unusually
small value of fsed, compared to the other planets’ best-fit value, while more
typical values of 0.1 and larger underestimate the observation. The CARMA
model albedo is somewhat lower, but is still within $2\sigma$ of the observed
value. If we wish to adapt our models to better match this planet’s high
dayside albedo, it likely would require increasing the predicted cloud
opacity, such as by increasing the porosity of the cloud particles (Samra et
al., 2020). Although increasing the atmospheric metallicity might also
increase the cloud opacity, published models for other planets indicate that
there is not a simple scaling between these two quantities (e.g. Morley et
al., 2013; Gao et al., 2018), and increasing the metallicity will also affect
the global thermal structure (e.g. Kataria et al., 2015).

Our model predictions for Kepler-17b and Kepler-8b are also in reasonable
agreement with the observed albedos. In both cases, an $f_{sed}=0.1$ Virga
model slightly overestimates the albedo while our CARMA model slightly
underestimates the albedo. This may indicate that moderately bright/cloudy
worlds have moderate sedimentation efficiencies. For Kepler-8b, both CARMA and
Virga models predict that the planet will have relatively uniform cloud
coverage in both latitude and longitude; we therefore do not need to consider
further spatial variations in cloud number density and particle size. For
Kepler-17b, the $f_{sed}=0.1$ Virga model predicts an albedo gradient across
the dayside atmosphere, but this gradient appears to be localized near the
equatorial (low latitude) region of the atmosphere. Our CARMA model predicts a
relatively uniform albedo across the two zones, but this may be biased by our
inability to resolve latitudinal gradients in the simplified two-zone model.

##  5 Conclusions

Optical secondary eclipse measurements made by _Kepler_ reveal a wide range of
geometric albedos for hot Jupiters with equilibrium temperatures between
1550-1700 K. We combine 3D general circulation models with both equilibrium
(Virga) and microphysical (CARMA) cloud models to explore whether 3D effects
can explain these observations. We find that the predicted albedos from our
Virga models are very sensitive to the assumed sedimentation efficiency
($f_{sed}$). We can compare these albedo predictions to results from our CARMA
model, which use mixing rates calculated from the GCM models to predict the
vertical extent and particle size distributions of the clouds. We find that
while the hemisphere-integrated CARMA albedos generally agree with the range
of albedos predicted by Virga, there is no single $f_{sed}$ value that
consistently matches the CARMA predictions.

When we compare these model predictions to the measured _Kepler_ albedos for
each of the six planets in our sample, we find that the albedos of K2-31b and
K2-107b are best matched by models that are either cloud-free or have very
deep compact cloud layers (large $f_{sed}$ values). Kepler-8b and Kepler-17b’s
optical albedos can be matched by moderately cloudy models ($f_{sed}$ greater
than 0.3). Both Virga and CARMA tend to under-predict the dayside albedos of
the two most reflective planets in our sample, HATS-11b and Kepler-7b, which
are best matched by Virga models with reflective Mg2SiO4 clouds extending to
very low pressures ($f_{sed}$ = 0.03); our CARMA model for HATS-11b predicts a
slightly lower albedo value than the brightest Virga model, while our CARMA
model for Kepler-7b predicts a slightly higher albedo value than Virga.
Although HATS-11b has relatively uniform cloud coverage across the dayside, it
is possible that other factors (such as a low particle porosity) might
increase the dayside cloud opacity beyond the values predicted by our models.

Our models predict that the observed albedo of Kepler-7b should be lower than
that of HATS-11b, in good agreement with the observations. Although a bright
reflective cloud layer forms in the westernmost region of the dayside
atmosphere, most dayside longitudes in Kepler-7b’s atmosphere are too warm for
Mg2SiO4 to condense, resulting in a lower hemisphere-averaged dayside albedo.
Empirical constraints on the horizontal extent of the western cloudy region
from phase curve observations indicate that it extends farther east than
predicted by our models, hinting that a more detailed study of the planet that
couples cloud microphysics and dynamics is required.

We conclude that the sample of optical albedos measured by _Kepler_ represents
a rich source of information for 3D cloud models, and that there is no single
explanation for the observed diversity of albedos for the planets considered
in this study. Future studies leveraging the large sample of transmission
spectra of hot Jupiters could provide complementary constraints on the typical
sedimentation efficiencies their atmospheres, while additional complementary
modeling studies exploring the coupled effects of atmospheric dynamics and
cloud microphysics, as well as an exploration of the micro-porosity of cloud
particles, would help to further illuminate the relative importance of these
processes in explaining the high albedos of the brightest planets in our
sample.

We thank the anonymous reviewer for an interesting report that greatly helped
improve the paper. P. Gao acknowledges support from NASA through the NASA
Hubble Fellowship grant HST-HF2-51456.001-A awarded by the Space Telescope
Science Institute, which is operated by the Association of Universities for
Research in Astronomy, Inc., for NASA, under contract NAS5-26555.

## References

  * Ackerman et al. (1995) Ackerman, A. S., Hobbs, P. V., & Toon, O. B. 1995, Journal of Atmospheric Sciences, 52, 1204 , doi: 10.1175/1520-0469(1995)052<1204:AMFPMT>2.0.CO;2
  * Ackerman & Marley (2001) Ackerman, A. S., & Marley, M. S. 2001, The Astrophysical Journal, 556, 872, doi: 10.1086/321540
  * Ackerman et al. (1993) Ackerman, A. S., Toon, O. B., & Hobbs, P. V. 1993, Science, 262, 226, doi: 10.1126/science.262.5131.226
  * Adams et al. (2019) Adams, D., Gao, P., de Pater, I., & Morley, C. V. 2019, The Astrophysical Journal, 874, 61, doi: 10.3847/1538-4357/ab074c
  * Adcrot et al. (2004) Adcrot, A., Campin, J.-M., Hill, C., & Marshall, J. 2004, Monthly Weather Review, 132, doi: 10.1175/MWR2823.1
  * Angerhausen et al. (2014) Angerhausen, D., Huber, K., Mandell, A., et al. 2014, Formation, Detection, and Characterization of Extrasolar Habitable Planets, Proceedings of the International Astronomical Union, IAU Symposium, 293, doi: 10.1017/S1743921313013343
  * Arakawa et al. (1977) Arakawa, E. T., Williams, M. W., & Inagaki, T. 1977, Journal of Applied Physics, 48, 3176, doi: 10.1063/1.324057
  * Barclay et al. (2012) Barclay, T., Huber, D., Rowe, J. F., et al. 2012, The Astrophysical Journal, 761, 53, doi: 10.1088/0004-637x/761/1/53
  * Barstow et al. (2014) Barstow, J. K., Aigrain, S., Irwin, P. G. J., et al. 2014, The Astrophysical Journal, 786, doi: 10.1088/0004-637X/786/2/154
  * Batalha et al. (2020) Batalha, N., caoimherooney11, & sagnickm. 2020, natashabatalha/virga: Initial Release, v0.0, Zenodo, doi: 10.5281/zenodo.3759888.  https://doi.org/10.5281/zenodo.3759888
  * Batalha et al. (2019) Batalha, N. E., Marley, M. S., Lewis, N. K., & Fortney, J. J. 2019, ApJ, 878, 70, doi: 10.3847/1538-4357/ab1b51
  * Bayliss et al. (2018) Bayliss, D., Hartman, J. D., Zhou, G., et al. 2018, The Astronomical Journal, 155, 119, doi: 10.3847/1538-3881/aaa8e6
  * Bell et al. (2017) Bell, T. J., Nikolov, N., Cowan, N. B., et al. 2017, The Astrophysical Journal, 847, L2, doi: 10.3847/2041-8213/aa876c
  * Bokeh Development Team (2014) Bokeh Development Team. 2014, Bokeh: Python library for interactive visualization
  * Bonomo, A. S. et al. (2017) Bonomo, A. S., Desidera, S., Benatti, S., et al. 2017, A&A, 602, A107, doi: 10.1051/0004-6361/201629882
  * Burrows et al. (2008) Burrows, A., Ibgui, L., & Hubeny, I. 2008, The Astrophysical Journal, 682, doi: 10.1086/589824
  * Chachan et al. (2019) Chachan, Y., Knutson, H. A., Gao, P., et al. 2019, The Astronomical Journal, 158, 244, doi: 10.3847/1538-3881/ab4e9a
  * Chachan et al. (2020) Chachan, Y., Jontof-Hutter, D., Knutson, H. A., et al. 2020, AJ, 160, 201, doi: 10.3847/1538-3881/abb23a
  * Colaprete et al. (1999) Colaprete, A., Toon, O. B., & Magalhães, J. A. 1999, Journal of Geophysical Research: Planets, 104, 9043, doi: https://doi.org/10.1029/1998JE900018
  * Coughlin & Lopez-Morales (2012) Coughlin, J., & Lopez-Morales, M. 2012, The Astrophysical Journal, 143, doi: 10.1088/0004-6256/143/2/39
  * Crossfield & Kreidberg (2017) Crossfield, I. J. M., & Kreidberg, L. 2017, The Astronomical Journal, 154, 261, doi: 10.3847/1538-3881/aa9279
  * Demory et al. (2011) Demory, B.-O., Seager, S., Madhusudhan, N., et al. 2011, The Astrophysical Journal, 735, doi: 10.1088/2041-8205/735/1/L12
  * Demory et al. (2013) Demory, B.-O., de Wit, J., Lewis, N., et al. 2013, The Astrophysical Journal Letters, 776, doi: 10.1088/2041-8205/776/2/L25
  * Desert et al. (2011) Desert, J.-M., Charbonneau, D., Fortney, J. J., et al. 2011, The Astrophysical Journal Supplement, 197, doi: 10.1088/0067-0049/197/1/11
  * Eigmuller et al. (2017) Eigmuller, P., Gandolfi, D., Persson, C. M., et al. 2017, The Astronomical Journal, 153, doi: 10.3847/1538-3881/aa5d0b
  * Esteves et al. (2013) Esteves, L. J., De Mooij, E. J. W., & Jayawardhana, R. 2013, The Astrophysical Journal, 772, doi: 10.1088/0004-637X/772/1/51
  * Esteves et al. (2015) —. 2015, The Astrophysical Journal, 804, doi: 10.1088/0004-637X/804/2/150
  * Evans et al. (2013) Evans, T. M., Pont, F., Sing, D. K., et al. 2013, The Astrophysical Journal Letters, 772, doi: 10.1088/2041-8205/772/2/L16
  * Fortney (2005) Fortney, J. J. 2005, Monthly Notices of the Royal Astronomical Society, 364, 649, doi: 10.1111/j.1365-2966.2005.09587.x
  * Fortney et al. (2006) Fortney, J. J., Cooper, C. S., Showman, A. P., Marley, M. S., & Freedman, R. S. 2006, The Astrophysical Journal, 652, 746, doi: 10.1086/508442
  * Fortney et al. (2011) Fortney, J. J., Demory, B.-O., Désert, J.-M., et al. 2011, The Astrophysical Journal Supplement Series, 197, 9, doi: 10.1088/0067-0049/197/1/9
  * Fraine et al. (2021) Fraine, J., Mayorga, L. C., Stevenson, K. B., et al. 2021, The Dark World: A Tale of WASP-43b in Reflected Light with HST WFC3/UVIS.  https://arxiv.org/abs/2103.16676
  * Fu et al. (2017) Fu, G., Deming, D., Knutson, H., et al. 2017, The Astrophysical Journal, 847, L22, doi: 10.3847/2041-8213/aa8e40
  * Gao et al. (2018) Gao, P., Marley, M. S., & Ackerman, A. S. 2018, The Astrophysical Journal, 855, 86, doi: 10.3847/1538-4357/aab0a1
  * Gao et al. (2017) Gao, P., Marley, M. S., Zahnle, K., Robinson, T. D., & Lewis, N. K. 2017, AJ, 153, 139, doi: 10.3847/1538-3881/aa5fab
  * Gao et al. (2014) Gao, P., Zhang, X., Crisp, D., Bardeen, C. G., & Yung, Y. L. 2014, Icarus, 231, 83, doi: https://doi.org/10.1016/j.icarus.2013.10.013
  * Gao et al. (2017) Gao, P., Fan, S., Wong, M. L., et al. 2017, Icarus, 287, 116 , doi: https://doi.org/10.1016/j.icarus.2016.09.030
  * Gao et al. (2020) Gao, P., Thorngren, D. P., Lee, G. K. H., et al. 2020, Nature Astronomy, 4, doi: 10.1038/s41550-020-1114-3
  * Goode et al. (2001) Goode, P. R., Qiu, J., Yurchyshyn, V., et al. 2001, Geophysical Research Letters, 28, 1671, doi: https://doi.org/10.1029/2000GL012580
  * Goody et al. (1989) Goody, R., West, R., Chen, L., & Crisp, D. 1989, Journal of Quantitative Spectroscopy and Radiative Transfer, 42, 539 , doi: https://doi.org/10.1016/0022-4073(89)90044-7
  * Grziwa et al. (2016) Grziwa, S., Gandolfi, D., Csizmadia, S., et al. 2016, The Astronomical Journal, 152, doi: 10.3847/0004-6256/152/5/132
  * Guillot & Showman (2002) Guillot, T., & Showman, A. P. 2002, Astronomy and Astrophysics, 385, doi: 10.1051/0004-6361:20011624
  * Helling (2019) Helling, C. 2019, Annual Review of Earth and Planetary Sciences, 47, doi: 10.1146/annurev-earth-053018-060401
  * Helling et al. (2008) Helling, C., Dehn, M., Woitke, P., & Hauschildt, P. H. 2008, The Astrophysical Journal, 675, L105, doi: 10.1086/533462
  * Helling et al. (2020) Helling, C., Kawashima, Y., Graham, V., et al. 2020, A&A, 641, A178, doi: 10.1051/0004-6361/202037633
  * Helling et al. (2019) Helling, C., Iro, N., Corrales, L., et al. 2019, Astronomy and Astrophysics, 631, doi: 10.1051/0004-6361/201935771
  * Heng & Demory (2013) Heng, K., & Demory, B.-O. 2013, The Astrophysical Journal, 777, 100, doi: 10.1088/0004-637x/777/2/100
  * Heng et al. (2021) Heng, K., Morris, B., & Kitzmann, D. 2021, arXiv, doi: arXiv:2103.02673
  * James et al. (1997) James, E., Toon, O., & Schubert, G. 1997, Icarus, 129, 147, doi: https://doi.org/10.1006/icar.1997.5763
  * Jensen & Toon (1994) Jensen, E. J., & Toon, O. B. 1994, Geophysical Research Letters, 21, 2019, doi: https://doi.org/10.1029/94GL01287
  * Kataria et al. (2015) Kataria, T., Showman, A. P., Fortney, J. J., et al. 2015, The Astrophysical Journal, 801, 86, doi: 10.1088/0004-637x/801/2/86
  * Kataria et al. (2013) Kataria, T., Showman, A. P., Lewis, N. K., et al. 2013, The Astrophysical Journal, 767, 76, doi: 10.1088/0004-637x/767/1/76
  * Kataria et al. (2016) Kataria, T., Sing, D. K., Lewis, N. K., et al. 2016, The Astrophysical Journal, 821, 9, doi: 10.3847/0004-637x/821/1/9
  * Kawashima & Ikoma (2019) Kawashima, Y., & Ikoma, M. 2019, The Astrophysical Journal, 877, 109, doi: 10.3847/1538-4357/ab1b1d
  * Kipping & Bakos (2011) Kipping, D., & Bakos, G. 2011, The Astrophysical Journal, 733, doi: 10.1088/0004-637X/733/1/336
  * Kirkpatrick (2005) Kirkpatrick, J. D. 2005, Annual Review of Astronomy and Astrophysics, 43, doi: 10.1146/annurev.astro.42.053102.134017
  * Kluyver et al. (2016) Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in ELPUB, 87–90
  * Lam et al. (2015) Lam, S. K., Pitrou, A., & Seibert, S. 2015, in Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC, LLVM ’15 (New York, NY, USA: ACM), 7:1–7:6.  http://doi.acm.org/10.1145/2833157.2833162
  * Lee et al. (2015) Lee, G., Helling, C., & Dobbs-Dixon, I. aand Juncher, D. 2015, Astronomy and Astrophysics, 580, doi: 10.1051/0004-6361/201525982
  * Lewis et al. (2014) Lewis, N. K., Showman, A. P., Fortney, J. J., Knutson, H. A., & Marley, M. S. 2014, The Astrophysical Journal, 795, 150, doi: 10.1088/0004-637x/795/2/150
  * Lines et al. (2019) Lines, S., Mayne, N. J., Manners, J., et al. 2019, Monthly Notices of the Royal Astronomical Society, 488, doi: 10.1093/mnras/stz1788
  * Lines et al. (2018) Lines, S., Mayne, N. J., Boutle, I. A., et al. 2018, Astronomy and Astrophysics, 615, doi: 10.1051/0004-6361/201732278
  * Livingston et al. (2018) Livingston, J. H., Crossfield, I. J. M., Petigura, E. A., et al. 2018, The Astronomical Journal, 156, 277, doi: 10.3847/1538-3881/aae778
  * Lodders (2003) Lodders, K. 2003, The Astrophysical Journal, 591, 1220, doi: 10.1086/375492
  * Lodders & Fegley (2002) Lodders, K., & Fegley, B. 2002, Icarus, 155, 393, doi: https://doi.org/10.1006/icar.2001.6740
  * MacDonald et al. (2018) MacDonald, R. J., Marley, M. S., Fortney, J. J., & Lewis, N. K. 2018, ApJ, 858, 69, doi: 10.3847/1538-4357/aabb05
  * Marley & Robinson (2015) Marley, M., & Robinson, T. 2015, Annual Review of Astronomy and Astrophysics, 53, 279, doi: 10.1146/annurev-astro-082214-122522
  * Marley et al. (2013) Marley, M. S., Ackerman, A. S., Cuzzi, J. N., & Kitzmann, D. 2013, Clouds and Hazes in Exoplanet Atmospheres, ed. S. J. Mackwell, A. A. Simon-Miller, J. W. Harder, & M. A. Bullock, 367
  * Marley et al. (1999) Marley, M. S., Gelino, C., Stephens, D., Lunine, J. I., & Freedman, R. 1999, The Astrophysical Journal, 513, 879, doi: 10.1086/306881
  * Marley et al. (2010) Marley, M. S., Saumon, D., & Goldblatt, C. 2010, The Astrophysical Journal Letters, 723, doi: 10.1088/2041-8205/723/1/L117
  * McGouldrick & Toon (2007) McGouldrick, K., & Toon, O. B. 2007, Icarus, 191, 1, doi: https://doi.org/10.1016/j.icarus.2007.04.007
  * McKinney (2010) McKinney, W. 2010, in Proceedings of the 9th Python in Science Conference, ed. S. van der Walt & J. Millman, 51 – 56
  * Mendonça et al. (2018) Mendonça, J. M., Malik, M., Demory, B.-O., & Heng, K. 2018, The Astronomical Journal, 155, 150, doi: 10.3847/1538-3881/aaaebc
  * Morley et al. (2013) Morley, C. V., Fortney, J. J., Kempton, E. M.-R., et al. 2013, The Astrophysical Journal, 775, 33, doi: 10.1088/0004-637x/775/1/33
  * Morley et al. (2012) Morley, C. V., Fortney, J. J., Marley, M. S., et al. 2012, The Astrophysical Journal, 756, 172, doi: 10.1088/0004-637x/756/2/172
  * Morley et al. (2015) —. 2015, The Astrophysical Journal, 815, 110, doi: 10.1088/0004-637x/815/2/110
  * Muñoz & Isaak (2015) Muñoz, A. G., & Isaak, K. G. 2015, Proceedings of the National Academy of Sciences, 112, 13461, doi: 10.1073/pnas.1509135112
  * Niraula et al. (2018) Niraula, P., Redfield, S., de Wit, J., et al. 2018, arXiv preprint arXiv:1812.09227
  * Oreshenko et al. (2016) Oreshenko, M., Heng, K., & Demory, B.-O. 2016, Monthly Notices of the Royal Astronomical Society, 457, doi: 10.1093/mnras/stw133
  * Parmentier et al. (2016) Parmentier, V., Fortney, J. J., Showman, A. P., Morley, C., & Marley, M. S. 2016, The Astrophysical Journal, 828, doi: 10.3847/0004-637X/828/1/22
  * Parmentier et al. (2021) Parmentier, V., Showman, A. P., & Fortney, J. J. 2021, MNRAS, 501, 78, doi: 10.1093/mnras/staa3418
  * Parmentier et al. (2013) Parmentier, V., Showman, A. P., & Lian, Y. 2013, A&A, 558, A91, doi: 10.1051/0004-6361/201321132
  * Parviainen et al. (2013) Parviainen, H., Deeg, H. J., & Belmonte, J. A. 2013, Astronomy & Astrophysics, 550, doi: 10.1051/0004-6361/201220081
  * Patzer et al. (2005) Patzer, A., Chang, C., Sedlmayr, E., & Sulzle, D. 2005, The European Physical Journal D - Atomic, Molecular, Optical and Plasma Physics, 32, doi: 10.1140/epjd/e2005-00026-8
  * Pérez & Granger (2007) Pérez, F., & Granger, B. E. 2007, Computing in Science & Engineering, 9
  * Powell et al. (2019) Powell, D., Louden, T., Kreidberg, L., et al. 2019, The Astrophysical Journal, 887, 170, doi: 10.3847/1538-4357/ab55d9
  * Powell et al. (2018) Powell, D., Zhang, X., Gao, P., & Parmentier, V. 2018, The Astrophysical Journal, 860, 18, doi: 10.3847/1538-4357/aac215
  * Pruppacher & Klett (1978) Pruppacher, H. R., & Klett, J. D. 1978, Homogeneous Nucleation (Dordrecht: Springer Netherlands), 162–182.  https://doi.org/10.1007/978-94-009-9905-3_7
  * Roman & Rauscher (2019) Roman, M., & Rauscher, E. 2019, The Astrophysical Journal, 872, doi: 10.3847/1538-4357/aafdb5
  * Roman et al. (2021) Roman, M. T., Kempton, E. M. R., Rauscher, E., et al. 2021, ApJ, 908, 101, doi: 10.3847/1538-4357/abd549
  * Rowe et al. (2008) Rowe, J., Matthews, J., Seager, S., et al. 2008, The Astrophysical Journal, 689, doi: 10.1086/591835
  * Samra et al. (2020) Samra, D., Helling, Ch., & Min, M. 2020, A&A, 639, A107, doi: 10.1051/0004-6361/202037553
  * Saumon & Marley (2008) Saumon, D., & Marley, M. S. 2008, The Astrophysical Journal, 689, 1327, doi: 10.1086/592734
  * Seager et al. (2000) Seager, S., Whitney, B. A., & Sasselov, D. D. 2000, The Astrophysical Journal, 540, doi: 10.1086/309292
  * Showman et al. (2009) Showman, A. P., Fortney, J. J., Lian, Y., et al. 2009, The Astrophysical Journal, 699, 564, doi: 10.1088/0004-637x/699/1/564
  * Showman & Kaspi (2013) Showman, A. P., & Kaspi, Y. 2013, The Astrophysical Journal, 776, 85, doi: 10.1088/0004-637x/776/2/85
  * Showman et al. (2015) Showman, A. P., Lewis, N. K., & Fortney, J. J. 2015, The Astrophysical Journal, 801, 95, doi: 10.1088/0004-637x/801/2/95
  * Showman et al. (2020) Showman, A. P., Tan, X., & Parmentier, V. 2020, Space Science Reviews, 216, doi: 10.1007/s11214-020-00758-8
  * Shporer & Hu (2015) Shporer, A., & Hu, R. 2015, The Astronomical Journal, 150, doi: 10.1088/0004-6256/150/4/112
  * Sing et al. (2016) Sing, D., Fortney, J., Nikolov, N., & et al. 2016, Nature, 529, 59, doi: https://doi.org/10.1038/nature16068
  * Skemer et al. (2016) Skemer, A. J., Morley, C. V., Zimmerman, N. T., et al. 2016, ApJ, 817, 166, doi: 10.3847/0004-637X/817/2/166
  * Smith (1998) Smith, M. D. 1998, Icarus, 132, 176, doi: https://doi.org/10.1006/icar.1997.5886
  * sqlite3 Development Team (2019) sqlite3 Development Team. 2019, DB-API 2.0 interface for SQLite databases
  * Sudarsky et al. (2003) Sudarsky, D., Burrows, A., & Hubeny, I. 2003, The Astrophysical Journal, 588, 1121, doi: 10.1086/374331
  * Toon et al. (1992) Toon, O., McKay, C., Griffith, C., & Turco, R. 1992, Icarus, 95, 24 , doi: https://doi.org/10.1016/0019-1035(92)90188-D
  * Visscher et al. (2010) Visscher, C., Lodders, K., & Fegley, B. 2010, The Astrophysical Journal, 716, 1060, doi: 10.1088/0004-637x/716/2/1060
  * von Paris et al. (2016) von Paris, P., Gratier, P., Borde, P., & Selsis, F. 2016, Astronomy & Astrophysics, 587, doi: 10.1051/0004-6361/201526297
  * Walt et al. (2011) Walt, S. v. d., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science & Engineering, 13, 22
  * Webber et al. (2015) Webber, M. W., Lewis, N. K., Marley, M., et al. 2015, The Astrophysical Journal, 804, doi: 10.1088/0004-637X/804/2/94
  * Wolf & Toon (2010) Wolf, E. T., & Toon, O. B. 2010, Science, 328, 1266, doi: 10.1126/science.1183260