Cloud Parameterizations and their Effect on Retrievals of Exoplanet Reflection Spectroscopy

Sagnick Mukherjee; Natasha E. Batalha; Mark S. Marley

doi:10.3847/1538-4357/abe53b

1. Introduction

Reflected light spectroscopy of exoplanets will be important in the next decade. Future space-based direct imaging missions will perform space-based optical coronagraphy, which will enable low-resolution (R ∼ 40–200) optical spectra of directly imaged cool giant planets and temperate terrestrial ones orbiting Sun-like stars. Although the Hubble Space Telescope (HST) or James Webb Space Telescope may be able to obtain transmission spectroscopy of a small number of favorable targets (e.g., HIP 41378 f; Dressing et al. 2020), reflected light offers the opportunity to directly infer the scattering properties of the atmosphere. A reflected light spectrum contains key information related to the planet's scattering and chemical properties. Specifically, unlike thermal and transmission spectroscopy, clouds and hazes determine the zeroth-order structure of the reflection spectra (Marley et al. 1999; Lupu et al. 2016; Gao et al. 2017; Nayak et al. 2017; MacDonald et al. 2018). Given the diversity of clouds and hazes expected in exoplanets (Morley et al. 2014; Gao et al. 2020), this presents new challenges for retrieving properties from reflected light spectra of exoplanets.

1.1. Previous Parameterizations Used in Reflected Light Retrievals

In order to prepare for this next decade of exoplanet spectroscopy there has been a growing body of literature aimed at determining best practices in retrieving properties from reflected light observations. We discuss each modeling framework below, focusing on the methodologies for parameterizing clouds. Given our focus on determining the required cloud complexity needed for retrievals, it is important to understand the parameterization method of each of these previous works.

1.1.1. Lupu et al. (2016) and Nayak et al. (2017)

Lupu et al. (2016) showed that the presence/absence of clouds and CH₄ can be inferred with high confidence from cool giants with reflection spectroscopy. The forward model used by Lupu et al. (2016) was based on the model initially developed by McKay et al. (1989), Marley & McKay (1999), and Marley et al. (1999), and later updated by Cahoy et al. (2010). They used CH₄ molecular opacity and the collision-induced opacities of CH₄, H₂, and He as gaseous opacities in their forward model because they worked only with planets where the reflection spectrum is expected to be dominated by CH₄. Lupu et al. (2016) used two simple retrieval models differing in cloud parameterization for retrieving on the reflected spectra, both of which contained wavelength-independent clouds. Their first model included a single semi-infinite cloud layer. Their second cloud model had a second cloud deck in addition to the semi-infinite cloud deck. In this case, the semi-infinite cloud deck at the bottom is forced to be optically thick and essentially acts as a reflective surface where the asymmetry parameter for this bottom deck is not retrieved. Using these two models, Lupu et al. (2016) retrieved on three validation test cases where the simulated data were produced using the retrieval model itself. They also tested their retrieval model on self-consistently modeled test cases of three planets. The performance of their retrieval models showed a decline on the "real" cases compared to the validation cases. For example, only lower limits on the CH₄ abundance could be obtained for two of the three real planet cases compared to constraints of CH₄ within factors of 10 of the true abundance for retrievals on validation planets. Lupu et al. (2016) hence note the need to (1) test various other parameterizations with varying levels of cloud modeling complexity on high signal-to-noise ratio (S/N) data of self-consistent models, and (2) identify an optimal set of cloud parameters that fully describe the system but minimize the number of free parameters. Using an identical modeling framework, Nayak et al. (2017) expanded on this work by exploring the effect of phase and radius uncertainty on the retrieval of atmospheric properties for cool giants with CH₄-dominated reflection spectra. Both of these studies concluded that an optical spectrum with a minimum S/N of 20 is required to retrieve accurate atmospheric properties such as molecular abundances and clouds from reflected light for cool giants.

1.1.2. Feng et al. (2018)

Feng et al. (2018) demonstrated the ability to ascertain the atmospheric composition of Earth-analog planets using higher-resolution spectroscopy (R = 70 or 140), as expected from mission concepts such as HabEx and LUVOIR. They considered only a single deck of H₂O cloud characterized by a cloud top pressure, cloud thickness, and a single value of optical depth. The asymmetry parameter and single scattering albedo in their model are fixed to constant values appropriate for H₂O clouds. Abundances of O₃, O₂, H₂O, and N₂ are retrieved, given the focus on Earth-like planet atmospheres. In addition to these parameters, they also retrieved a parameter for the reflective surface and a parameter to describe the patchiness of clouds, f_c. Although their methodology pertains to terrestrial planets, the patchy cloud concept will likely be relevant for gas giants as well.

In order to incorporate patchy cloud coverage, they create two models for each case: one with 100% cloud coverage and a second that is cloud-free. The albedo spectra of each of the runs are combined by weighting the first cloudy spectrum with f_c and the latter cloudless spectrum with (1 – f_c). They found that at relatively higher spectral resolutions of R ∼ 140 and S/N ∼ 20 the parameters of interest could be constrained. At R ∼ 70 and S/N ∼ 20 the presence of clouds and molecules could be detected, and with low resolution (R ∼ 50) combined with S/N ∼ 20 one could achieve just weak detections of clouds and molecules. This retrieval model only accounts for water clouds. Terrestrial planet atmospheres might have photochemical hazes, which would complicate inferring scattering properties with water-only single cloud deck models. Lastly, Feng et al. (2018) kept their single scattering albedo and asymmetry parameter fixed. This motivates additional work to determine how retrieving on these parameters for single or multiple cloud decks affects the retrieved solutions for cool giants.

1.1.3. Hu (2019) and Damiano & Hu (2020)

A separate model, ExoREL, was developed by Hu (2019) for modeling reflected spectra of cool giants, and was implemented in a retrieval framework by Damiano & Hu (2020). This model considered the effect of H₂O and NH₃ condensation on the vapor-phase mixing ratio of these molecules.

Damiano & Hu (2020) conducted the retrieval analysis on three test planets, including 47 UMa b, a focus of this analysis as well. Similar to Lupu et al. (2016) and Nayak et al. (2017), Damiano & Hu (2020) retrieved the cloud bottom pressure, cloud thickness, a mixing ratio for well-mixed CH₄, and the gravity. But they assumed only H₂O and NH₃ condensation and focused on retrieving the volume mixing ratio (VMR) of H₂O and NH₃ below the cloud bottom. They retrieve a condensation ratio, which is then used to construct the depleted VMR of CH₄ and NH₃ above the cloud deck. Damiano & Hu (2020) retrieved on the albedo spectra of three test planets synthesized using the forward model of Hu (2019).

1.1.4. Irwin et al. (2008) and Barstow et al. (2014)

NEMESIS is a well-vetted code for retrieving properties of solar system planets (Irwin et al. 2008) and exoplanets (Barstow et al. 2014). Most recently, it was used to demonstrate the ability to retrieve cloud scattering properties in thermal emission (Taylor et al. 2020). With regard to reflected light, Barstow et al. (2014) used NEMESIS to determine the atmospheric parameters from the hot-Jupiter HD 189733 b, observed using HST/Space Telescope Imaging Spectrograph (STIS) by Evans et al. (2013) at very low spectral resolution. Instead of a Bayesian retrieval analysis, they perform a chi-square analysis on a grid of 980 spectra with a fixed set of cloud base pressures, particle sizes, and optical depths at 0.25 μm—and a fixed set of chemistry parameters for the VMR of Na. Similar to previous analyses of cool giants where the cloud species is assumed, Barstow et al. (2014) assume that the clouds are composed of MgSiO₃ or MnS because these are the relevant condensation species for HD 189733 b. Then the scattering properties for the cloud particles were calculated using Mie theory with a double-peaked Henyey–Greenstein formulation of the phase function. Other important atmospheric parameters such as the temperature–pressure profile and VMRs of CO, CH₄, H₂O, and CO₂ were fixed at the best-fit values derived from previously obtained emission spectra. They found that the data were consistent with a large number of cloudy cases, as well as many cloudless cases. Therefore, there was lot of degeneracy in their cloud parameter space. Given the data quality and resolution of STIS, this is consistent with the findings of Lupu et al. (2016), Nayak et al. (2017), Feng et al. (2018), and Damiano & Hu (2020), who reported the need for an S/N of 20 for proper characterization of cloud properties.

NEMESIS has also been used for retrievals of scattering properties of solar system planets from reflection spectra (e.g., Irwin et al. 2015, 2016). Unlike in the study of exoplanets, the cloud species, location, and thickness for solar system planets are generally known quantities and can be treated as fixed parameters. Moreover, the data quality used in these solar system studies is generally far superior in quality to the focus of this and other previous studies discussed thus far. Therefore in studies such as Irwin et al. (2015), the wavelength-dependent imaginary part of the refractive indices and the parameters of the particle size distribution for the each condensing species can be directly retrieved from the data.

1.2. This Analysis

Previous work has highlighted the need for retrieval studies to be conducted on self-consistent models without the assumption of pure water clouds. Therefore, in this work, we build upon the approach of Lupu et al. (2016) and others by retrieving the atmospheric properties of cool giant planets from spectra produced with a cloud model. We test four parameterizations to retrieve the cloud properties, each with an increasing level of complexity. Ultimately, we try to compare the accuracy and precision of the retrieved chemistry and clouds of the atmospheres to the original input. We focus on three radiative properties of the clouds—optical depth per layer, asymmetry parameter, and the single scattering albedo. By retrieving directly on the cloud radiative properties, we avoid any assumptions about the condensing species in the atmosphere. For our simulated data, we model the reflected spectra for three cool giants that are priority targets for the Nancy Grace Roman Space Telescope. These priority cool giant targets have high contrast ratios compared to their host stars in reflected light due to the optimal combination of their size, separation from host star, and effective temperature.

Table 1 summarizes the spectral resolutions in the optical wavelength range expected from multiple future space-based direct imaging missions (The LUVOIR Team 2019; Gaudi et al. 2020). We consider both the lowest spectral resolution expected from the Nancy Grace Roman Space Telescope and higher spectral resolutions expected from mission concepts such as HabEx and LUVOIR while producing mock observations of the reflected spectra of these three priority target exoplanets. We retrieve on these mock observation spectra in order to test and compare various methodologies for parameterizing atmospheres when retrieving properties. This exercise helps to inform the complexity of atmospheric parameterization requisite for the next decade of reflected light studies. We also account for the uncertainty in the gravity of the planets while performing the retrievals. In doing this, we aim to address the following:

1.
Does the choice of cloud parameterization affect the cloud properties and molecular mixing ratios retrieved from the reflected spectra of cool giant planets?
2.
Does the performance of our retrieval model change from one planet to another (i.e., different effective temperatures)?
3.
Does the constraint on the retrieved gravity depend on the cloud parameterization?
4.
How does data quality (S/N and R) limit ability to retrieve molecular abundances and cloud properties?

Table 1. Expected Spectral Resolutions from Future Space-based Direct Imaging Missions

Future Mission	Wavelength Range (μm)	Spectral Resolution
Roman Space Telescope CGI Spectroscopy	0.675–0.785	47–75^a
Roman Space Telescope CGI Imaging	0.5–0.8^a	⋯
LUVOIR A & B (ECLIPS)	0.515–1.03	140^b
HabEx Coronagraph	0.45–1.00	140^c
HabEx Starshade Instrument	0.45–0.975	140^c

Notes.

^a https://roman.ipac.caltech.edu/sims/Param_db.html.^bThe LUVOIR Team (2019).^cGaudi et al. (2020).

Download table as: ASCII Typeset image

We describe the atmospheric modeling using PICASO and Virga in Section 2. We briefly explain the retrieval method in Section 3. We describe the results of our analysis in Section 4 and finally discuss our main results in Section 5 and summarize them in Section 6.

2. Modeling Reflected Spectra

We use the parameter space of effective temperature versus gravity to select a representative target population for the analysis. We calculate the equilibrium temperature (T_eq) and gravity (g) of 23 direct imaging planet targets, most of which have radial velocity detections (Fischer et al. 2002; Butler et al. 2006; Hatzes et al. 2006; Howard & Fulton 2016). We use the planetary orbital and stellar parameters to calculate the equilibrium temperature of these planets assuming a zero albedo. We also consider an additional internal temperature of 100 K, similar to that of Jupiter (Fortney et al. 2007), to get the effective temperature (T_eff). We use the M sin i and planet radius for calculating the gravity of these planets. We calculate the planet radius using the empirical mass–radius relationship for cool giants from Thorngren et al. (2019). The T_eff versus g parameter space for these planets is shown in Figure 1. We use eps Eri b, 47 UMa b, and HD 62509 b as our target planets (shown in Figure 1) to explore retrievals on reflected light across three different temperature regimes for directly imaged cool giants. We emphasize that the aim is not to produce highly self-consistent models of these planets. Instead, we aim to explore a range in temperature that enables a diversity in cloud formation and chemistry scenarios. Additionally, our targets have similar gravity estimates. This allows us to isolate the effect that varying cloud and chemistry scenarios has on retrieving atmospheric parameters. In a future analysis we will explore the effect of completely unconstrained gravity.

**Figure 1.** Effective temperature T_eff vs. gravity g parameter space for 23 cool giant targets for space-based direct imaging missions. The blue, orange, and green points denote the target planets used in this work: eps Eri b, 47 UMa b, and HD 62509 b, respectively. The upward arrows point to the effective temperature if a 100 K internal temperature is considered. Main point—We choose three cool giant target planets that span a wide range of effective temperatures (135–533 K) in order to explore diversity in atmospheric cloud properties.
Download figure:
Standard image High-resolution image

2.1. Modeling the Planet Atmospheres

We use PICASO (Batalha et al. 2019), which has heritage from McKay et al. (1989), Marley & McKay (1999), and Cahoy et al. (2010), to model the reflected light spectra of our target planets. PICASO is an open-source radiative transfer code capable of calculating transmitted, reflected, and/or thermal spectra of planets and brown dwarfs. PICASO requires the temperature–pressure (T(P)) profile, cloud structure and atmospheric chemistry as inputs for the radiative transfer calculation. Here, we discuss modeling each of these inputs for our reflected spectra simulation in Sections 2.1.1–2.1.3 and discuss the basics of PICASO in Section 2.1.4.

2.1.1. The Temperature–Pressure Profile

We divide the planet atmosphere into 61 plane-parallel pressure layers where the pressure rises logarithmically from 10⁻⁶ to 10³ bar. We model the temperature–pressure profile of the planet atmospheres using the empirical parameterization described in MacDonald et al. (2018). This empirical T(P) profile parameterization is dependent on T_eff, g, and the metallicity of the planet [M/H]. The best-fit values of the coefficients in the parameterization have been determined by fitting the empirical profile to a large number of T(P) profiles for cool giants produced self-consistently using the methodology of Fortney et al. (2008). The T(P) profile parameterization is described by

$\begin{eqnarray}&&{T}^{4}(P)={T}_{0}^{4}+{T}_{{\rm{d}}{\rm{e}}{\rm{e}}{\rm{p}}}^{4}\left(\displaystyle \frac{P}{1000\,{\rm{b}}{\rm{a}}{\rm{r}}}\right)\end{eqnarray} \tag{ 1 }$

where both T₀ and T_deep are functions dependent on T_eff, g, and [M/H] (see MacDonald et al. 2018). We assume a Jupiter metallicity of 3× solar metallicity (Wong et al. 2004) for all the three planets. The parameterized T(P) profiles for the three planets are shown in Figure 2. Although the self-consistent T(P) profiles will show more structure than the parameterized profiles, as seen in MacDonald et al. (2018), this will not affect the results of the analysis given the insensitivity of reflected light spectroscopy to temperature.

**Figure 2.** Temperature–pressure profiles of eps Eri b (blue), 47 UMa b (orange), and HD 62509 b (green) generated using the empirical parameterization of MacDonald et al. (2018). Condensation curves of the molecular condensates considered in the cloud calculation of `Virga` are shown by the colored dashed lines. Main point—Our target planets have different T(P) profiles, which causes condensation of different species at different pressures, resulting in varying cloud decks and optical properties.
Download figure:
Standard image High-resolution image

2.1.2. Atmospheric Chemistry

Chemical equilibrium abundances are interpolated from those computed on a grid of (P, T) points as calculated using a modified version of the NASA CEA Gibbs minimization code (see Gordon & McBride 1994). The chemistry grid is available for download from Marley et al. (2018) and described fully in M. S. Marley et al. (2021, in preparation). For these cool giants the most important gaseous absorbers in the optical are methane, ammonia, and water. The grid accounts for depletion of each chemical species above the point of condensation. The VMR profiles of these three species for each of our three target planets are shown in Figure 3.

**Figure 3.** Volume mixing ratio as a function of pressure of CH₄, NH₃, and H₂O for the three planets. Main point—The three planets have very different chemical structure due to different T(P) profiles.
Download figure:
Standard image High-resolution image

2.1.3. Clouds with `Virga`

We calculate the cloud profiles for each of the targets using Virga. Virga follows the treatment of condensation in atmospheres by Ackerman & Marley (2001). For each condensate species and at each atmospheric layer, the vapor pressure in excess of the saturation vapor pressure is allowed to condense. The condensation curves of all the molecular species considered by our cloud model are shown in Figure 2 in dashed lines. The condensation curves follow those in Morley et al. (2012) and Gao et al. (2020), and are available online.⁴

The full cloud profiles from Virga are shown Figure 4. The case of the coolest planet, eps Eri b, has H₂O and NH₃ condensation forming two separate cloud decks. The warmer case, 47 UMa b, lacks NH₃ clouds but is still dominated by H₂O clouds. The hottest case, HD 62509 b, is dominated by Na₂S clouds at depth (P > 1 bar), and lacks condensation from H₂O or NH₃. As shown in Figure 4, the three cases explored here probe three different cloud condensation regimes.

The vertical structure of the cloud is determined by the balance between the vertical turbulent mixing of condensates and vapor and the sedimentation of condensates described by the equation

$\begin{eqnarray}&&-{K}_{z}\displaystyle \frac{\partial {q}_{t}}{\partial z}-{f}_{\mathrm{sed}}{w}_{* }{q}_{c}=0\end{eqnarray} \tag{ 2 }$

where q_c is the condensate mole fraction, q_v is the vapor mole fraction, and q_t = q_c + q_v (Ackerman & Marley 2001). The first term represents the vertical turbulent mixing of the condensate and vapor, where K_z is the vertical eddy diffusion coefficient. The second term represents the sedimentation caused by the condensates, where w_* is the convective velocity scale. f_sed is a dimensionless ratio of the sedimentation velocity to w_*. The cloud structure is solved by the balance of the two competing processes for each of the condensing species.

The sedimentation parameter f_sed and the vertical eddy diffusion coefficient K_z are the two inputs that critically determine the vertical extents and particle size distributions of the clouds. Overall, low values of f_sed < 1 produce thick cloud layers with smaller particles, and higher values, f_sed > 1, produce thinner cloud decks with larger particles. We set f_sed = 3, motivated by the higher values that have been successful in modeling cool giant clouds for Jupiter-like planets, in contrast to lower values of ∼0.1, which have typically been used for hot-Jupiters (Webber et al. 2015).

The vertical eddy diffusion coefficient (K_z) can strongly affect the cloud properties. Overall, higher K_z values lead to the formation of larger cloud particles. We calculate the vertical eddy diffusion coefficient K_z using (Gierasch & Conrath 1985)

$\begin{eqnarray}&&{K}_{z}=\displaystyle \frac{H}{3}{\left(\displaystyle \frac{L}{H}\right)}^{4/3}{\left(\displaystyle \frac{{RF}}{\mu {\rho }_{a}{c}_{p}}\right)}^{1/3}\end{eqnarray} \tag{ 3 }$

where H is the atmospheric scale height, L is the turbulent mixing length, R is the universal gas constant, F is the thermal flux of the atmosphere (assumed to be $\sigma {T}_{\mathrm{eff}}^{4}$ ), μ is the atmospheric molecular weight, ρ_a is the atmospheric density, and c_p is the atmospheric specific heat at constant pressure.

This formulation is based on the assumption that the vertical eddy diffusion coefficient for the vapor and the condensate of the cloud model is the same as derived for heat in free convection conditions (Gierasch & Conrath 1985). This method assumes that convection occurs all the way to the top of the atmosphere, which of course is not the case in reality. However, this has been used to baseline the cloud model for Jupiter in Ackerman & Marley (2001) and hence is applicable for this study.

In solving Equation (2), we compute an effective particle radius per layer per species, and assume a log-normal distribution of particles with a geometric standard deviation of 2 about that radius. The Mie scattering calculations are then computed with PyMieScatt (Sumlin et al. 2018) over this distribution. This allows Virga to produce the altitude- and wavelength-dependent optical depth per layer τ(P, λ), single scattering albedo ω(P, λ), and asymmetry parameter g(P, λ) of the clouds in the atmosphere. The single scattering albedo describes the wavelength-dependent reflectivity of the cloud particles. Higher ω leads to higher reflectivity. The asymmetry parameter captures the forward scattering/backscattering probability of the scattering of light from the cloud particles.

The final cloud optical profiles of each planet are shown in Figure 5. The optical properties can trace back to the exact cloud species. For example, the ∼1 bar cloud deck of eps Eri b corresponds to the highest region of single scattering and therefore can be reasonably identified as an H₂O cloud. Ultimately, it is the information in these profiles that we aim to recover.

2.1.4. Reflected Spectra with `PICASO`

With the T(P) profile, chemical structure, and the cloud optical properties as inputs, we adopt the one-dimensional version of PICASO to calculate the reflected light spectra for the three planets (Batalha et al. 2019). PICASO uses the two-stream radiative transfer methodology of Toon et al. (1989). We add contributions from both molecular and collision-induced absorption (CIA) opacities. Although we only focus on retrieving H₂O, NH₃, and CH₄, PICASO includes the molecular opacity from H₂O (Barber et al. 2006; Tennyson & Yurchenko 2018), CH₄ (Yurchenko et al. 2013; Yurchenko & Tennyson 2014), NH₃ (Yurchenko et al. 2011), CO (Li et al. 2015), PH₃ (Sousa-Silva et al. 2014), H₂S (Azzam et al. 2016), CO₂ (Huang et al. 2014), Na and K (Ryabchikova et al. 2015), and others not applicable to these temperatures (e.g., TiO, VO). Among the CIA PICASO includes opacity from H₂–H₂ (Abel et al. 2011), H₂–He, H₂–N₂, H₂–H, H₂–CH₄, H–electron bound–free, H–electron free–free, and H₂–electron interactions. The resultant opacity calculations are available on Zenodo (Batalha et al. 2020b).

In order to accurately capture asymmetrical backscattering caused by Rayleigh scattering, we use the two-term Henyey–Greenstein (TTHG) phase function combined with Rayleigh phase function formalism (TTHG_Ray in PICASO) for the direct scattering component. The one-term Henyey–Greenstein (P_OTHG) phase function is defined as

$\begin{eqnarray}&&{P}_{\mathrm{OTHG}}(\cos {\rm{\Theta }})=\displaystyle \frac{1}{2}\displaystyle \frac{1-{g}^{2}}{\sqrt{1+{g}^{2}-2g\cos {\rm{\Theta }}}}.\end{eqnarray} \tag{ 4 }$

The TTHG phase function capturing both forward scattering, g_f, and backscattering, g_b, is then defined as

$\begin{eqnarray}\begin{array}{rcl}{P}_{\mathrm{TTHG}}(\cos {\rm{\Theta }}) & = & {{fP}}_{\mathrm{OTHG}}(\cos {\rm{\Theta }},{g}_{f})\\ & & +(1-f){P}_{\mathrm{OTHG}}(\cos {\rm{\Theta }},{g}_{b})\end{array}\end{eqnarray} \tag{ 5 }$

where g_f = $\bar{g}$ , g_b = – $\bar{g}$ /2, and f, the fraction of forward scattering to backscattering, is $f=1-{g}_{b}^{2}$ . The first term corresponds to the forward scattering phase function weighted by f, while the second term is the backscattering phase function weighted by (1 − f). Our adoption of g_b is arbitrary. However, it has been previously adopted in studies of exoplanet reflected light (Cahoy et al. 2010; Feng et al. 2018) due to the lack of a priori information. $\bar{g}$ is the cloud asymmetry parameter weighted by the cloud fractional opacity. It is calculated using the cloud asymmetry parameter, single scattering albedo, cloud opacity, and Rayleigh opacity:

$\begin{eqnarray}&&\bar{g}=\left(\displaystyle \frac{\omega {\tau }_{\mathrm{cld}}}{\omega {\tau }_{\mathrm{cld}}+{\tau }_{\mathrm{Ray}}}\right)g.\end{eqnarray} \tag{ 6 }$

The TTHG and Rayleigh scattering phase functions (P_Ray) are then combined with a weighted addition to get the final phase function, TTHG_Ray. The weight factors for the TTHG and Rayleigh phase functions are τ_cld/τ_scat and τ_Ray/τ_scat, respectively. P_TTHGRay is then

$\begin{eqnarray}&&{P}_{{\rm{TTHGRay}}}(\cos {\rm{\Theta }})=\displaystyle \frac{{\tau }_{{\rm{cld}}}}{{\tau }_{{\rm{scat}}}}{P}_{{\rm{TTHG}}}(\cos {\rm{\Theta }})+\displaystyle \frac{{\tau }_{{\rm{Ray}}}}{{\tau }_{{\rm{scat}}}}{P}_{{\rm{Ray}}}(\cos {\rm{\Theta }}).\end{eqnarray} \tag{ 7 }$

The multiple scattering phase function in PICASO is calculated by expanding the HG function to second order (N=2) and forcing the second-order moment such that it reproduces Rayleigh scattering when it dominates the total scattering opacity (Batalha et al. 2019). We also include the effect of Raman scattering by using the formalism of Pollack et al. (1986). Their methodology for Raman scattering results in redshift of the photons, which reduces the overall reflectively toward the blue (Batalha et al. 2019). The formulation of Pollack et al. (1986) does not model high resolution of solar emission features seen in reflected light spectra of gas giants (Oklopčić et al. 2016). However, these features require far too high a resolution (R ≫ 100) for the next decade of direct imaging observations (Oklopčić et al. 2016).

In order to understand the interplay between Rayleigh, molecular, and cloud scattering, PICASO computes the "photon attenuation," which denotes the pressure level where the two-way optical depth from each component reaches τ = 1. Figure 6 shows the photon attenuation for the three planets. The "flatness" of the reflected spectra of eps Eri b comes from the dominance of the optical properties of water clouds, which are highly reflective and non-wavelength-dependent at the respective particle radii from 0.3 to 1 μm. On the other hand, the HD 62509 b spectrum is dominated by Rayleigh opacity short of 0.5 μm, and by molecular absorption between 0.5 and 1 μm. Because the cloud deck is much lower (in altitude) than the molecular opacity source, the molecular absorption dominates, causing the planet to have the lowest albedo among the three planets. Molecular, cloud, and Rayleigh opacities contribute significantly to the reflected light of 47 UMa b.

For the case of eps Eri b and 47 UMa b, which both have contributions from the cloud, the different behavior in the spectra is a result of different cloud optical properties (optical depth, single scattering, and asymmetry profiles). As seen from Figures 5 and 4, eps Eri b has two cloud decks (an NH₃ cloud deck above a very optically thick H₂O cloud deck), whereas 47 UMa b only has an H₂O cloud deck at a pressure similar to the higher NH₃ cloud deck in eps Eri b. The H₂O cloud deck in eps Eri b has a relatively high optical depth (∼200) compared to the H₂O cloud deck in 47 UMa b. 47 UMa b also has a second cloud deck, which appears much deeper in the atmosphere and is relatively optically thick compared to the deeper ZnS cloud deck appearing in eps Eri b. The optical depth differences, in addition to the differing asymmetry and the single scattering albedo of the H₂O cloud deck of eps Eri b, lead to a flat cloud-dominated albedo spectrum for eps Eri b compared to that for 47 UMa b.

Finally, the simulated reflected light spectra for the three planets are shown in the lower panel of Figure 6. They show (1) a case dominated by bright water clouds, (2) a case dominated by clouds, molecular opacity, and Rayleigh scattering, and lastly (3) a case dominated by Rayleigh and molecular opacity. These three spectra are exemplary cases to be used in the retrieval analysis because they test the parameterizations under three different scattering regimes.

2.2. The Retrieval Setup

Similar to previous works discussed, a retrieval requires parameterizations to be made in order to best capture the behavior of the physical model, described in Section 2.1. For the retrieval, we replace the chemistry, cloud, and T(P) profile with parameterizations that can be used in a Bayesian framework. We use an isothermal T(P) profile for our forward model with the temperature fixed at the effective temperature of the planet. This is different from Lupu et al. (2016), where the T(P) profile for retrievals was kept fixed to the profile used for modeling the simulated data. Unlike thermal and transmission spectroscopy, reflected light is only sensitive to the temperature through its contribution to the line shapes of the molecular opacity and through the scale height of the atmosphere. The opacity is not strongly temperature-dependent in the parameter space probed by our three targets (Karkoschka 1994; Karkoschka & Tomasko 2011).

We initially assume that the atmosphere is well mixed and start by retrieving a single value for the mixing ratios of three molecules—CH₄, H₂O, and NH₃—that dominate the opacity sources for this cooler class of planets (Burrows et al. 1997; Madhusudhan et al. 2016). The rest of the atmospheric composition, other than CH₄, NH₃, and H₂O, is assumed to be composed of H₂ and He. The H₂/He fraction is taken to be f = H₂/He = 0.837/0.163 (Lodders 2019). Therefore the He and H₂ mixing ratio is given by

$\begin{eqnarray}&&\mathrm{He}=\displaystyle \frac{1-{\mathrm{CH}}_{4}-{{\rm{H}}}_{2}{\rm{O}}-{\mathrm{NH}}_{3}}{f+1},\end{eqnarray} \tag{ 8 }$

$\begin{eqnarray}&&{{\rm{H}}}_{2}=f\times {\rm{He}}\end{eqnarray} \tag{ 9 }$

where the species' name represents its VMRs (v/v). Our initial assumption of well-mixed profiles in the retrieval is different from that of Damiano & Hu (2020), who use two free parameters for H₂O and NH₃ each to describe their depleted mixing ratios. Our chemical profiles incorporate depletion caused by condensation, as is clearly evident in Figure 3. eps Eri b shows a depletion of NH₃ and H₂O due to condensation of both species whereas 47 UMa b shows a depletion only in the H₂O mixing ratio because NH₃ condensation is absent from 47 UMa b. Therefore, our initial assumption of well-mixed atmospheres purposely tests whether or not additional complexity is necessary. Later, we relax this assumption and discuss the results of a retrieved depleted profile in Section 4.1.2.

We parameterize τ(P), ω(P), and g(P) in four ways with different levels of complexity. We neglect any wavelength dependence in all the three cloud optical properties for our retrieval model. The validity of this assumption, given that Virga calculates wavelength-dependent cloud optical properties, is addressed in Section 5. A schematic diagram for the cloud parameterizations is shown in Figure 7. In what follows, we describe each of the four parameterizations (9–15 free parameters in total), and the associated priors used in the retrieval analysis.

**Figure 7.** Schematic diagram of the four cloud parameterizations of structure and optics used in our forward model for the target planets. The left column shows the optical depth per layer, τ; the middle column is for the asymmetry parameter, g, per layer; and the right column shows the single scattering albedo, ω, per layer. The top row shows the parameterizations for Case 1—the box cloud model, the second row is for Case 2—the single cloud profile model, the third row shows Case 3—the double cloud profile model, and the last row shows Case 4—the double cloud profile model with two-valued asymmetry and the single scattering albedo model. Details of the model parameters for each case can be found in Sections 2.2.1 through 2.2.4.
Download figure:
Standard image High-resolution image

2.2.1. The Box Cloud Model (Case 1)

The optical depth profile (τ(P)) of this cloud parameterization is similar to the cloud deck used for retrieval by Feng et al. (2018). Unlike Feng et al. (2018), the asymmetry parameter and the single scattering albedo for this parameterization are also free parameters. We model the cloud structure in this case using five parameters. The other four parameters are the CH₄, NH₃, and H₂O mixing ratios and the gravity of the planet. The cloud structure for this case is parameterized according to the following equations:

$\begin{eqnarray}\tau (P,\lambda )=\left\{\begin{array}{ll}0 & \mathrm{if}\ P\gt {P}_{0}\ \mathrm{or}\ P\lt {P}_{0}-{dP},\\ {\tau }_{\max } & \mathrm{otherwise}\end{array}\right.\end{eqnarray} \tag{ 10 }$

$\begin{eqnarray}\omega (P,\lambda )=\left\{\begin{array}{ll}0 & \mathrm{if}\ P\gt {P}_{0}\ \mathrm{or}\ P\lt {P}_{0}-{dP},\\ \omega & \mathrm{otherwise}\end{array}\right.\end{eqnarray} \tag{ 11 }$

$\begin{eqnarray}g(P,\lambda )=\left\{\begin{array}{ll}0 & \mathrm{if}\ P\gt {P}_{0}\ \mathrm{or}\ P\lt {P}_{0}-{dP},\\ g & \mathrm{otherwise}\end{array}\right.\end{eqnarray} \tag{ 12 }$

where P₀, ${\tau }_{\max }$ , ω, g, and dP are the five parameters of the model. The nine parameters of this model and the priors used in retrievals are summarized in Table 2.

Table 2. Parameters for the Box Cloud Model

Parameter	Description	Range	Iteration type
g	Asymmetry parameter	0–1	Linear
ω	Single scattering albedo	0–1	Linear
τ	Optical depth per layer	0.1–30	Linear
P₀	Cloud base pressure level	10⁻⁶–10² bar	Log scale
dP	Cloud deck thickness	0 to P₀ – 10⁻⁶ bar	Log scale
CH₄	CH₄ mixing ratio	−6 to 0	Log scale
NH₃	NH₃ mixing ratio	−6 to 0	Log scale
H₂O	H₂O mixing ratio	−6 to 0	Log scale
g	Gravity	25–65 m s⁻²	Linear

Download table as: ASCII Typeset image

2.2.2. Single Cloud Profile Model (Case 2)

For this case, we model the clouds with an altitude-dependent optical depth profile. The asymmetry parameter and the single scattering albedo are forced to be zero beneath the base of the cloud deck and can take a value between zero and one above the base of the cloud deck up to the top of the atmosphere. The cloud structure here is modeled as

$\begin{eqnarray}\tau (P,\lambda )=\left\{\begin{array}{ll}0 & \mathrm{if}\ P\gt {P}_{0},\\ {\tau }_{\max }{e}^{-a{\left(\mathrm{ln}P-\mathrm{ln}{P}_{0}\right)}^{2}} & \mathrm{otherwise}\end{array}\right.\end{eqnarray} \tag{ 13 }$

$\begin{eqnarray}\omega (P,\lambda )=\left\{\begin{array}{ll}0 & \mathrm{if}\ P\gt {P}_{0},\\ \omega & \mathrm{otherwise}\end{array}\right.\end{eqnarray} \tag{ 14 }$

$\begin{eqnarray}g(P,\lambda )=\left\{\begin{array}{ll}0 & \mathrm{if}\,P\gt {P}_{0},\\ g & \mathrm{otherwise}\end{array}\right.\end{eqnarray} \tag{ 15 }$

where P₀, ${\tau }_{\max }$ , a, ω, and g are free parameters. Hence, our forward model involves nine free parameters consisting of five cloud parameters and the mixing ratios and gravity, similar to Case 1. The number of parameters is the same as in Case 1 and they are described in Table 3.

Table 3. Parameters for the Single Cloud Profile Model

Parameter	Description	Range	Iteration type
g	Asymmetry parameter	0–1	Linear
ω	Single scattering albedo	0–1	Linear
τ	Optical depth per layer	0.1–30	Linear
P₀	Cloud base pressure	10⁻⁶–10² bar	Log scale
a	Cloud deck scale height	10⁻⁴ to 2	Log scale
CH₄	CH₄ mixing ratio	−6 to 0	Log scale
NH₃	NH₃ mixing ratio	−6 to 0	Log scale
H₂O	H₂O mixing ratio	−6 to 0	Log scale
g	Gravity	25–65 m s⁻²	Linear

Download table as: ASCII Typeset image

2.2.3. Double Cloud Profile Model (Case 3)

This model is similar to Case 2 except a second cloud deck is allowed to form here. Additionally, the asymmetry and single scattering profiles are similar to Case 2, where they take a value between zero and one above the base of the deepest cloud deck. The 12 free parameters are described in Table 4. The following equations describes the parameterizations for this model:

$\begin{eqnarray}\begin{array}{c}\tau (P,\lambda )=\left\{\begin{array}{ll}0 & {\rm{i}}{\rm{f}}\,P\gt {P}_{1},\\ {\tau }_{max1}{e}^{-{a}_{1}{\left({\rm{l}}{\rm{n}}P-{\rm{l}}{\rm{n}}{P}_{1}\right)}^{2}} & {\rm{i}}{\rm{f}}\,P\lt {P}_{1},\\ {\tau }_{max2}{e}^{-{a}_{2}{\left({\rm{l}}{\rm{n}}P-{\rm{l}}{\rm{n}}{P}_{2}\right)}^{2}} & {\rm{i}}{\rm{f}}\,P\lt {P}_{2},\end{array}\right.\end{array}\end{eqnarray} \tag{ 16 }$

$\begin{eqnarray}\omega (P,\lambda )=\left\{\begin{array}{ll}0 & \mathrm{if}\ P\gt {P}_{1},\\ \omega & \mathrm{otherwise}\end{array}\right.\end{eqnarray} \tag{ 17 }$

$\begin{eqnarray}g(P,\lambda )=\left\{\begin{array}{ll}0 & \mathrm{if}\ P\gt {P}_{1},\\ g & \mathrm{otherwise}\end{array}\right.\end{eqnarray} \tag{ 18 }$

Table 4. Parameters for the Double Cloud Deck Profile Model

Parameter	Description	Range	Iteration Type
g	Asymmetry parameter	0–1	Linear
ω	Single scattering albedo	0–1	Linear
τ₁	Optical depth per layer of lower cloud deck	0.1–30	Linear
τ₂	Optical depth per layer of upper cloud deck	0.1–30	Linear
P₁	Cloud base pressure of lower cloud deck	10⁻⁶—10² bar	Log scale
P₂	Cloud base pressure of upper cloud deck	P₁–10⁻⁶ bar	Log scale
a₁	Cloud deck scale height of lower cloud deck	10⁻⁴ to 2	Log scale
a₂	Cloud deck scale height of upper cloud deck	10⁻⁴ to 2	Log scale
CH₄	CH₄ mixing ratio	−6 to 0	Log scale
NH₃	NH₃ mixing ratio	−6 to 0	Log scale
H₂O	H₂O mixing ratio	−6 to 0	Log scale
g	Gravity	25–65 m s⁻²	Linear

Download table as: ASCII Typeset image

2.2.4. Double Cloud Profile Model with Two-valued g₀ and w₀ (Case 4)

This 15-parameter model has the same parameterization of optical depth per layer as Case 3. The major difference in this case is that the asymmetry and single scattering are each allowed to have two values. The asymmetry and single scattering parameterizations are

$\begin{eqnarray}\omega (P,\lambda )=\left\{\begin{array}{ll}{\omega }_{2} & \mathrm{if}\ P\gt {P}_{2}\ \mathrm{and}\ P\lt {P}_{2}-{dP},\\ {\omega }_{1} & \mathrm{otherwise}\end{array}\right.\end{eqnarray} \tag{ 19 }$

$\begin{eqnarray}\begin{array}{c}g(P,\lambda )=\left\{\begin{array}{ll}{g}_{2} & {\rm{i}}{\rm{f}}\,P\gt {P}_{2}\,{\rm{a}}{\rm{n}}{\rm{d}}\,P\lt {P}_{2}-dP,\\ {g}_{1} & {\rm{o}}{\rm{t}}{\rm{h}}{\rm{e}}{\rm{r}}{\rm{w}}{\rm{i}}{\rm{s}}{\rm{e}}.\end{array}\right.\end{array}\end{eqnarray} \tag{ 20 }$

Table 5 describes the 15 parameters for Case 4.

Table 5. Parameters for the Double Cloud Deck Profile Model with Two-valued g and w

Parameter	Description	Range	Iteration type
g₁	Asymmetry parameter	0–1	Linear
g₂	Asymmetry parameter	0–1	Linear
ω₁	Single scattering albedo	0–1	Linear
ω₂	Single scattering albedo	0–1	Linear
τ₁	Optical depth per layer of lower cloud deck	0.1–30	Linear
τ₂	Optical depth per layer of upper cloud deck	0.1–30	Linear
P₁	Cloud base pressure level of lower cloud deck	10⁻⁶–10² bar	Log scale
P₂	Cloud base pressure level of upper cloud deck	P₁ – 10⁻⁶ bar	Log scale
a₁	Cloud deck scale height of lower cloud deck	10⁻⁴ to 2	Log scale
a₂	Cloud deck scale height of upper cloud deck	10⁻⁴ to 2	Log scale
dP	Thickness of g₀/w₀ deck	0 to P₂ – 10⁻⁶ bar	Log scale
CH₄	CH₄ mixing ratio	−6 to 0	Log scale
NH₃	NH₃ mixing ratio	−6 to 0	Log scale
H₂O	H₂O mixing ratio	−6 to 0	Log scale
g	Gravity	25–65 m s⁻²	Linear

Download table as: ASCII Typeset image

3. Dynamic Nested Sampling

We use the Dynamic Nested Sampling package Dynesty ⁵ (Speagle 2020) for our retrievals. We choose Dynamic Nested Sampling to explore non-Gaussian posteriors of the parameters. The Nested Sampling method (Skilling 2006) can efficiently and accurately determine both the evidence and the posteriors of the problem simultaneously, unlike traditional Markov Chain Monte Carlo, which prioritizes the estimation of the posterior. Briefly, the basic process includes (1) drawing a large number (∼50 × number of free parameters) of live points from the priors of the parameter space provided, and then (2) iteratively replacing the live point with the least likelihood with a new live point having a greater likelihood than the replaced point. At each iteration, the replaced points become "dead points." The evidence can then be estimated with a set of N dead points by summing over the product of their likelihoods and prior volumes (Skilling 2006). This process continues until a user-defined stopping criterion is met.

We wrap the retrieval model described in Section 2.2 in the Dynesty module and retrieve on the simulated observational spectra. In each retrieval, we assign 50 live points per free parameter, as recommended (Speagle 2020). We use the the multi-ellipsoidal decomposition method because of its ability to efficiently capture complex, multi-modal posteriors (Speagle 2020). We use ( ${\rm{\Delta }}\mathrm{ln}({Z}_{i})$ ) defined as

$\begin{eqnarray}&&{\rm{\Delta }}\mathrm{ln}({Z}_{i})=\mathrm{ln}({Z}_{i}+{\rm{\Delta }}{Z}_{i})-\mathrm{ln}({Z}_{i})\end{eqnarray} \tag{ 21 }$

as our stopping criterion. Here, Z_i and ΔZ_i are the current and remaining evidence estimates, respectively. The remaining evidence can be approximated by the product of the highest likelihood among remaining live points and the prior volume of the last dead point. The retrieval is stopped when this ${\rm{\Delta }}\mathrm{ln}({Z}_{i})$ is smaller than N/1000, where N is the number of live points. This criterion has been optimized for evidence and posterior estimation by Speagle (2020). We do not use any limitation on the maximum number of iterations for the retrievals.

We use the evidences estimated from the nested sampling for calculating the Bayes factor. The Bayes factor allows us to quantitatively compare each of our models. This can directly inform us whether or not one model is favored over another. The Bayes factor of model M₀ over model M₁ for a data set D is

$\begin{eqnarray}&&{B}_{01}=\displaystyle \frac{p(D| {M}_{0})}{p(D| {M}_{1})}\end{eqnarray} \tag{ 22 }$

where p(D∣M₀) and p(D∣M₁) are the evidence of model 0 and of model 1 with the data set D, respectively. Pairs of models with $\mathrm{ln}(B)$ less than 2.1 are said to indicate that model 0 is weakly favored over model 1 at best, with a confidence of less than 2σ (Trotta 2008). If $\mathrm{ln}(B)$ is greater than 5, model 1 can be strongly ruled out over model 0 with a confidence of 3.6σ–5σ (Trotta 2008). We present the retrieval results obtained using Dynesty in the following sections.

4. Results

We first bin the model spectra for eps Eri b, 47 UMa b, and HD 62509 b to a constant resolution (R) of 40. We also fix the S/N of the spectra to 20 at a wavelength of 0.35 μm. This initial choice of S/N is motivated by earlier studies (Lupu et al. 2016; Nayak et al. 2017; Feng et al. 2018; Hu 2019) establishing this to be the minimum S/N required for accurate retrievals of atmospheric properties. The values for R and S/N are chosen to mimic the likely best possible data quality from space-based direct imaging and spectroscopy missions in the near future. We also explore the effect of degrading S/N to 5. We retrieve on the albedo spectra for each planet using our four parameterizations described in Section 2.2. Here we present the results of our analysis on a planet-by-planet basis.

4.1. 47 UMa b

4.1.1. Comparison of Retrieval Parameterizations

For the first case, we retrieve atmospheric parameters using all four cases described in Section 2.2 on 47 UMa b. We estimate the effective temperature of 47 UMa b to be ∼217 K. Figure 8 shows the median retrieved solution along with 1σ and 2σ confidence intervals. The residuals of the retrieved median spectra from the simulated data are also shown for each case in Figure 8.

**Figure 8.** The top panel shows the model albedo spectrum for 47 UMa b. The second, third, fourth, and fifth panels show the comparison of retrieved spectra by each retrieval model from Case 1 to Case 4 with the simulated observed spectrum of 47 UMa b shown as black diamonds. The red line shows the median albedo spectrum of the retrieved models while the dark and light blue shaded regions represent the 1σ and 2σ uncertainties, respectively. Main point—Cases 3 and 4 provide much better fits to the spectra than Cases 1 and 2.
Download figure:
Standard image High-resolution image

Cases 1 and 2 show greater residuals on the blue side of the spectra than Cases 3 and 4. The large residuals toward the blue for the retrievals in Cases 1 and 2 foreshadow a potential overestimation of the molecular abundances. However, comparing the performances of Case 3 and Case 4 with just the residuals and the median spectra is not quantitatively informative. Hence, we use the evidence estimates from the nested sampling calculations in order to determine how strongly one model could be ruled out or compared to another model. Specifically, we calculate the Bayes factor described in Section 3 for each pair of models. The heat map is shown in Figure 9. As suggested by Figure 8, we see that both Cases 3 and 4 are favored over Cases 1 and 2. Of the two nine-parameter models, Case 2 is favored very weakly over Case 1. The 15-parameter model of Case 4 is moderately favored over the 12-parameter model of Case 3. Moving forward, we evaluate each of the retrieval models by comparing their retrievals of various atmospheric properties with the input properties used to simulate the mock spectra.

**Figure 9.** Heat map of the natural logarithm of the Bayes factor of "Model 0" on the y-axis over "Model 1" on the x-axis. A higher Bayes factor, as is seen for Cases 3 and 4 compared to Cases 1 and 2, signifies that Cases 3 and 4 can be favored over Cases 1 and 2. Between Cases 1 and 2, Case 2 very weakly rejects Case 1. Between Case 3 and 4, Case 4 moderately rejects Case 3. Main point—Cases 3 and 4 outperform Cases 1 and 2 in retrievals. Case 4 (highest complexity) is moderately favored over Case 3.
Download figure:
Standard image High-resolution image

Figure 10 presents the comparison of the retrieved molecular mixing ratios compared with the input molecular mixing ratio profiles, which was used to generate the simulated observed spectra for 47 UMa b. In this first case, a single value for the mixing ratio was retrieved. The full posteriors of this single retrieved value are depicted in light blue, while the altitude-dependent profile is shown in dashed red.

H₂O is the most dominant molecular opacity source for 47 UMa b, among the three molecules. The H₂O opacity causes broadband molecular absorption features starting from 0.5–1 μm. Constraints on the H₂O mixing ratio for Cases 1 and 2 are precise but not accurate. Constraints on H₂O mixing ratio for Cases 3 and Case 4 are relatively accurate, but the precision appears to be obscured by the altitude dependence. The posterior of the H₂O mixing ratios for Cases 3 and 4 is most highly peaked at the region of highest water mixing ratio of the "true" state of the atmosphere (deep within it). However, the posterior contains an additional tail that extends to the depleted value. This might indicate that there is sensitivity to the depleted abundances above the cloud deck. In Section 4.1.2 we report the results of adding in an additional parameter to capture the depleted water profile.

CH₄ is retrieved by both Case 3 and Case 4 within 2σ of the "true" CH₄ profile, though the magnitude of the 2σ value is relatively large (an order of magnitude in abundance) and the posterior is not Gaussian. For this planet, this result is not surprising. CH₄ contributes to some of the absorption features beyond 0.7 μm, but unlike eps Eri b the CH₄ opacity contribution is smaller than H₂O at all wavelengths. While Case 3 retrieves CH₄ well within the 2σ limit of the input CH₄ mixing ratio profile, Case 4 better constrains CH₄ within 1σ. Like H₂O, CH₄ is overestimated by Cases 1 and 2 with sharply peaked posteriors.

In all cases, S/N = 20 does not allow for the precise retrieval of the abundance of NH₃ because, although it is relatively abundant, its opacity contribution to the spectrum is negligible until 0.8 μm. Cases 1 and 2 fail to constrain NH₃ at all, and return nearly uniform posteriors (the prior). Case 3 also returns a posterior whose 1σ constraint spans ∼2 orders of magnitude. The retrieved posterior distribution for NH₃ with Case 4 shows a bias toward higher NH₃ abundance, and does not yield a Gaussian posterior, which points to some degeneracy. With the CH₄, H₂O, and NH₃ retrievals for 47 UMa b described above, it is clear that Cases 1 and 2 are not robust models for retrieving abundances while Cases 3 and 4 retrieve the mixing ratios of all the three molecules within 2σ of the true profiles.

Cloud retrievals for each case are shown in Figure 11. The dashed red lines show the 0.3–1 μm averaged input profiles produced by Virga. The left column shows the retrievals of optical depth per layer for the four cases. The middle and right columns show the retrieval of the asymmetry and single scattering albedo, respectively. The median profiles (in blue) along with 1σ and 2σ bounds (in dark and light brown) are shown for those cases where altitude dependence was considered. For the rest, the posteriors for the single parameter are shown by the blue shaded curves.

The profiles of optical depth per layer show that both the single-deck parameterizations—Cases 1 and 2—trace only the upper cloud deck. Case 1 retrieves the position of the top H₂O cloud deck more precisely than Case 2, whereas Case 2 retrieves the position of the top cloud deck more accurately than Case 1. Case 3 retrieves the position of the bottom cloud deck and optical depth profile accurately but fails to retrieve accurate or precise parameters for the top cloud deck. Similarly, Case 4 gets the bottom deck but places the top cloud deck at much lower pressures (<10⁻⁴ bar), missing the "true" cloud deck pressure by several (∼4) orders of magnitude.

The asymmetry parameters retrieved by Cases 1 and 2 are highly degenerate (the posterior is multi-modal). Case 3 retrieval of the profile of the asymmetry parameter is relatively precise but it overestimates its value. Case 4 retrieves an inaccurate asymmetry parameter profile with a low (∼0.2) asymmetry parameter across most of the atmosphere but a large value (∼0.9) at the position of the retrieved top cloud deck.

High values (∼1) of the single scattering albedo are retrieved by Cases 1, 2, and 3. These three retrieval models hence prefer reflective cloud particles, similar to H₂O. This is not the case for Case 4, where the retrieved bottom deck is relatively unreflective with a single scattering albedo of 0.4, whereas the top deck is highly reflective with a single scattering albedo close to 1.

None of the four retrieval models succeeds in retrieving the cloud structure for 47 UMa b completely. Some aspects of the position of the cloud layers are accurately captured by Cases 2, 3, and 4. For example, Cases 3 and 4 retrieve the position and the optical depth profile of the deeper bottom cloud deck accurately, but incorrectly retrieve either the position or optical depth of the top layer. We discuss this further in Section 5.

4.1.2. 47 UMa b Case 3 with a Depleted Water Profile

The posterior distribution of the well-mixed H₂O mixing ratio retrieved using Case 3 on 47 UMa b shows a bias toward lower abundances of H₂O, while the retrieved median value coincides with the H₂O mixing ratio below the cloud deck in the "true" atmosphere. This bias, shown in Figure 10, might be indicative of model sensitivity with respect to the depletion above the top H₂O cloud deck. This sensitivity was similarly highlighted in the work of Damiano & Hu (2020). Here, we modify the well-mixed H₂O retrieval model in Case 3 to accommodate a depleted water mixing ratio above the top cloud deck. This requires two additional parameters for Case 3 since we now retrieve two H₂O mixing ratios—one below (H₂O¹) and one above (H₂O²) the cloud deck. In the retrieval model, we define the cloud top pressure, P_top, as the location where the optical depth of the top cloud deck is 10^z times less than the peak cloud optical depth, where z is the second additional parameter. The water mixing ratio is modeled to be a straight line with negative slope in the log(pressure) versus log(H₂O) space between the lower and upper cloud deck pressures.

Figure 12 shows the spectral fit of the retrieved results along with the cloud and molecular abundance retrievals. We find that this additional complexity within Case 3 does not improve the retrieved cloud optical depth structure significantly for 47 UMa b compared to the well-mixed Case 3 retrieval. Without the depleted water profile, Case 3 is unable to constrain the top cloud deck pressure within the lower bounds of the cloud deck pressure parameter space but with the depleted water profile, the additional complexity in this modified Case 3 model helps to put 1σ bounds on the top cloud pressure within the lowest pressure bound of the retrieved atmospheric model. The retrieval of the deeper cloud deck remains similarly accurate/precise in both models. This parameterization also broadens the constraints on the asymmetry parameter toward lower values but the overestimation of the asymmetry parameter as seen for Case 3 in Figure 11 still persists. Retrieval of the single scattering albedo remains similar in both cases.

This model captures the depleted H₂O abundance above the cloud deck within 1σ and can also retrieve the deep H₂O abundance within 1σ. This modified model does not show any improvement in the retrieval of CH₄ over the original Case 3 model. NH₃ is better constrained by the modified Case 3 model than Case 3 (unmodified), where NH₃ is not at all constrained.

Despite the apparent improvement in water abundance, a Bayes factor analysis suggests that Case 3 without water depletion is weakly preferred over Case 3 with water depletion for 47 UMa b. The difference in Bayes factor is so small that without prior knowledge of the true solution, it would be difficult to discern which scenario was correct. That is, the complexity that arises from the addition of the depleted profile parameter is not strongly favored.

4.1.3. Retrieval at Higher Spectral Resolution (R = 140)

Future mission concepts such as LUVOIR and HabEx are expected to achieve a spectral resolution of ∼140 in optical wavelengths as shown in Table 1. We retrieve on an R ∼ 140 and S/N ∼ 20 spectrum of 47 UMa b with the double cloud deck model, Case 3, to investigate how a higher spectral resolution from these missions can improve upon the constraints on the molecular abundances and cloud properties obtained from data with lower spectral resolution expected from the Roman Space Telescope.

Figure 13 shows corner plots of the retrieved Case 3 parameters from the 47 UMa b reflection spectrum with spectral resolution of R ∼ 40 and 140 (left and right panels, respectively). Here, we opt to show the full set of marginalized posterior distributions of each parameter, in order to showcase correlations between the various retrieved parameters. Overall, a higher-resolution spectrum leads to tighter constraints on multiple parameters, relative to the R ∼ 40 spectrum with Case 3.

Specifically, the bias in the H₂O retrieval seen in the lower-resolution retrieval with Case 3 is no longer present in the posterior retrieved from the the higher-resolution spectrum. The 2σ constraints on both H₂O and CH₄ have tightened by a factor of ∼2 relative to the posteriors shown in Figure 13. CH₄ retrieval has improved from being within 2σ of the "true" CH₄ abundance in the low-resolution case to being within 1σ at high resolution. The NH₃ posteriors remain unconstrained at both spectral resolutions—though this is as expected given its opacity contribution at these wavelengths. Lastly, although the abundances improve significantly, retrievals of the cloud parameters remain very similar (within 1σ) to that obtained from the lower-resolution spectra. This further motivates the utility of low-resolution spectroscopy for obtaining cloud properties.

4.2. eps Eri b

Out of our sample of three planets, the spectroscopic model of eps Eri b is the most similar to Jupiter, and it is our coolest target. The reflected spectrum is dominated by cloud opacity, as seen in Figure 6. Therefore, it has a reflection spectrum of zero slope with a major CH₄ feature at 0.73 and 0.9 μm along with a minor NH₃ feature. The retrieval results are shown in Figure 14.

**Figure 14.** Comparison of the retrieved atmospheric properties of eps Eri b. The left column shows the retrieval results using Case 2 parameterization while the right column shows retrieval results with Case 3 parameterizations. The top row shows the comparison of the observed spectrum in black diamonds with the median (red), 1σ (dark blue), and 2σ (light blue) regions of the retrieved spectra with each parameterization. The middle row shows the retrieved mixing ratios of CH₄, NH₃, and H₂O from left to right in both columns with the true input profiles shown in red. The last row shows the retrieved cloud properties—optical depth per layer, asymmetry parameter, and single scattering albedo from left to right in both columns with the true input profiles shown in red. Main point—Case 2 performs slightly better than Case 3, and is preferred according to its Bayes factor.
Download figure:
Standard image High-resolution image

Overall, the accuracy and precision of the retrieved cloud and mixing ratio profiles of the atmosphere were found to be very similar for Cases 2 and 3, which negates the addition of more complexity in going to Case 4. Moreover, the Bayes factor calculation suggests that Case 2 is weakly favored over Case 3.

Despite the multiple cloud decks made up of different condensate species, Case 3 retrieves a single deck. The maximum optical depth of the highest pressure water cloud located at 1 bar is τ > 10. The NH₃ cloud directly above the water deck only reaches a maximum optical depth of τ > 0.5. Therefore, when retrieving the optical depth profile, the model favors a single larger cloud deck, which spans both the H₂O and the NH₃ clouds. Despite H₂O and NH₃ having different optical properties, the asymmetry parameter and single scattering albedo structure for the planet are retrieved accurately and precisely within the 1σ bounds of the input profile by both cases. Here, we define accuracy as retrieving the true value in the region of maximum cloud opacity.

Both Cases 2 and 3 overestimate the CH₄ abundance for eps Eri b despite the presence of two prominent CH₄ features in the eps Eri b albedo spectrum. This overestimation of CH₄ was also noted in Lupu et al. (2016) for the case of HD 99492 c, which also contained two CH₄ absorption features similar to the eps Eri b spectrum in our case. NH₃ is neither precise nor accurate, which is intuitive given the opacity contribution to the spectrum. H₂O remains unconstrained for eps Eri b because the "true" H₂O mixing ratio depletes at high pressures (>1 bar) due to the formation of water clouds. Therefore, there is a limited H₂O contribution to the albedo spectrum for eps Eri b.

We performed additional retrieval tests for eps Eri b in order to determine the source of the overestimated CH₄. The tests we performed included (1) fixing the gravity to the true value, (2) increasing the resolution of the simulated spectrum to R = 120, and (3) removing water as a free parameter. Fixing gravity and removal of H₂O did not increase the accuracy or precision of the CH₄ abundance. Increasing the resolution to R = 120 increased the accuracy of the retrieved CH₄ by two orders of magnitude, which is a clear indication that higher resolution is needed to accurately retrieve CH₄. At R = 120, four data points sample the undersaturated CH₄ feature at 0.7 μm, while eight data points sample the saturated feature at 0.9 μm. This is in contrast to R = 40, where only one data point samples the undersaturated feature at 0.7 μm and two sample the feature at 0.9 μm.

4.3. HD 62509 b

HD 62509 b is our hottest target planet with an estimated effective temperature of ∼533 K. The clouds in this planet are formed much deeper in the atmosphere than in the other two cooler targets in our consideration, as is evident in Figure 2. As a result, most of the spectrum is dominated by molecular opacity and Rayleigh scattering as shown in Figure 6. Additionally, the planet is dim in reflected light compared to the other two planets. Similar to eps Eri b, the Bayes factor calculation suggests that Case 2 is weakly preferred over Case 3 for HD 62509 b. Therefore, the accuracy and precision of both Case 2 and Case 3 in retrieving cloud properties and the molecular abundances are very similar for HD 62509 b.

H₂O mixing ratio is retrieved accurately within 2σ by Cases 2 and 3 for HD 62509 b as shown in Figure 15. This is because, like 47 UMa b, H₂O is the most significant gaseous opacity source in the atmosphere of HD 62509 b. Both Cases 2 and 3 fail to constrain the CH₄ and NH₃ abundances for HD 62509 b since the spectrum is dominated by H₂O opacity for HD 62509 b with negligible opacity contribution from CH₄ and NH₃.

For HD 62509 b, the model of Cases 2 and 3 traces the region where the lower cloud deck achieves optical depth of ∼10. The retrieval is able to constrain the cloud optical depth structure despite (1) the cloud being very deep in the atmosphere and (2) the relatively small contribution of cloud opacity to the reflected spectrum compared to our two other target planets. Both Cases 2 and 3 underestimate the single scattering albedo for HD 62509 b. This is likely because the cloud deck is at such depth in the atmosphere that molecular absorption dominates the total opacity. The retrieval models hence become less sensitive to the reflectivity of the cloud particles. The posterior of the asymmetry parameter nearly spans the entire prior region. This too is because of the lack of sensitivity of the model to the cloud scattering properties due to the depth of the cloud in the atmosphere.

5. Discussion

5.1. Required Cloud Complexity

Of the first two cases with nine parameters, Case 1 (box model) could never accurately/precisely retrieve atmospheric or scattering properties. Therefore, we advocate against box models for future retrieval work, even in simplified studies.

For Cases 2 and 3 (one or two cloud decks), the results were less clearly defined. The reflected light spectrum of eps Eri b is dominated by cloud opacity, while the hotter HD 62509 b spectrum is dominated by molecular and Rayleigh opacity. In both cases we find that a single-deck model (e.g., Case 2) can sufficiently capture the atmospheric properties and produce similar results to the double cloud deck model (e.g., Case 3). For the case of 47 UMa b, though, double deck models (Cases 3 and 4) are significantly favored over single-deck models (Cases 1 and 2). Naively, this is contrary to the intuition that might be gained from looking at Figure 4, which shows a strong NH₃ and CH₄ double cloud deck for eps Eri b.

The key factor that dictated whether an additional cloud deck was necessary was the pressure location of the cloud deck with the highest optical depth (τ ∼ 10)—relative to the region of maximum contribution of molecular opacity. The eps Eri b NH₃ cloud deck achieved τ ∼ 0.5 at 0.5 bar. The H₂O deck achieved τ ∼ 10 at 1 bar. The bulk of the molecular opacity (see Figure 6) was mostly contained below this optically thick H₂O deck (e.g., cloud layer–cloud layer–molecular opacity). This is in contrast to 47 UMa b, where the upper cloud deck achieved τ ∼ 1 at 0.1 bar while the lower deck did not achieve τ ∼ 10 until roughly 30 bar. The molecular opacity sat between these two regions (e.g., cloud–molecular–cloud). Therefore, the former case of eps Eri b could be accurately modeled with one larger cloud deck, whereas the latter case of 47 UMa b required two separate scattering regions, one below and one above the region of highest molecular opacity. Of course we will not know a priori where these cloud decks exist. However, this result can inform the interpretation of future reflected light results. Retrievals that strongly favor single-deck models may not ultimately reflect the true state.

Lastly, Case 4—two parameters to describe asymmetry and single scattering albedo—was only moderately favored over Case 3, according to Bayes factor analysis. However, Case 4 failed to retrieve an accurate and precise single scattering albedo and the asymmetry parameter profile for 47 UMa b. Additionally, the accuracy of the molecular abundances retrieved with Case 4 was also similar to Case 3 (within 1σ). Even though the use of Case 4 was not strongly motivated in this work, it could be useful for cases not explored here. Specifically, the use of Case 4 would be suitable for atmospheres that have stronger vertical variation in asymmetry, or single scattering albedo. For example, a case with a water cloud deck below high-altitude photochemically produced hazes (e.g., Gao et al. 2017) might require at least two parameters in each of the scattering properties.

5.1.1. Ability to Constrain Gravity

Unlike for most transiting planets, there generally will only be approximate constraints on gravity for directly imaged planets in reflected light. While radial velocity and a sufficient number of images will constrain $\sin i$ , planet radii will still be uncertain. Therefore it is worthwhile to determine (1) whether or not gravity can be accurately retrieved from reflected light spectroscopy alone, and (2) whether or not an imprecise gravity affects the ability to retrieve accurate atmospheric properties.

In order to determine the robustness of our results with respect to imprecise gravity measurements, we allow the gravity of our planets to vary by ±50% of the assumed mass. With this level of uncertainty we can then explore whether or not the knowledge of the mass can be improved with the observed 0.3–1 μm reflected light spectroscopy. Figure 16 shows the posteriors for the retrieved gravity for the five parameterizations explored for 47 UMa b in Section 4.1.1.

**Figure 16.** Comparison of the retrieved gravity posteriors for the four cases and modified Case 3 on 47 UMa b. The black dashed line depicts the true gravity of the planet. The green shaded curve depicts the posterior for Case 1 while the blue shaded curve depicts the posterior for Case 2. The posterior for Case 3 is depicted by the red shaded curve, modified Case 3 is depicted by the black shaded curve, and that for Case 4 is depicted with the blue shaded curve. Main point—No parameterization can accurately improve the gravity beyond ±50% of the a priori value.
Download figure:
Standard image High-resolution image

Cases 1 and 2 retrieve gravity posteriors that clearly favor higher gravity values whereas Case 3 (no water depletion), Case 3 (with water depletion), and Case 4 show posteriors peaking toward lower values of gravity (beyond 2σ) than the true input gravity shown by the dotted line in Figure 16. From this we draw two conclusions: (1) none of the parameterizations here could reliably retrieve gravity with the spectral resolution and S/N of reflected spectra used in this work; (2) even with ±50% inaccuracy in gravity cloud structures the abundances of molecules can be inferred from the reflected spectra of cool giants with proper choice of cloud parameterizations. Further analysis, beyond the scope of this analysis, would need to be performed to determine whether these conclusions were (1) robust against gravity constraints that were larger than ±50%, and (2) robust against unconstrained phase angle (see further discussion on phase in Section 5.3).

5.1.2. Effect of Assumed Signal-to-noise Ratio

Throughout the analysis we fixed S/N. We determined that for a planet such as 47 UMa b the retrieved results were highly dependent on the complexity of the retrieval model used and the overall parameterization. In order to determine the robustness of this result with respect to the assumed S/N, we degrade the S/N to see whether this complexity dependence still holds for a simulated spectrum of 47 UMa b with a lower S/N of 5.

Figure 17 shows the retrieval results on a spectrum with S/N = 5 for Case 2 and Case 3 parameterizations. At lower S/N, retrievals produced with Case 1 and Case 2 parameterizations are able to fit the observed spectra. This is an intuitive result because the extra absorption features, which were seen at S/N = 20, are now buried within the systematic error bars of the simulated spectra. Similar to previous results of Lupu et al. (2016) and Feng et al. (2018), we find that at such low S/N, none of the cases results in precise or accurate constraints on molecular abundances. According to the Bayes factor, Case 2 rules out Case 3 very weakly with this quality of data.

Despite not being able to constrain molecular abundances directly, we can make inferences as to where the cloud deck is relative to the molecular and Rayleigh opacity levels based on the retrieved photon attenuation map. We demonstrate this with the retrieved photon attenuation map for Case 2 shown in Figure 18. Comparing this retrieved photon attenuation map with that shown in Figure 6 shows that the cloud optical depth level can be retrieved within 2σ of the "true" opacity levels with the Case 2 parameterization. The Rayleigh opacity levels and the gas opacity levels are also retrieved within 1σ of their respective "true" optical depth. This estimate of the cloud base pressure level from an albedo spectra with S/N = 5 can roughly and indirectly inform the temperature–pressure structure of the atmosphere. This is because the location of the cloud deck is predicated on the region where the temperature becomes cool enough to condense a particular species. Therefore, by combining the expected equilibrium temperature of the planet and the retrieved cloud deck, zeroth-order inferences can be made about the potential temperature regime of the atmosphere.

**Figure 18.** Map of photon attenuation depth corresponding to an optical depth of 0.5 retrieved using Case 2 from the S/N = 5 spectrum of 47 UMa b. The attenuation pressure levels are divided into Rayleigh, gas, and cloud opacity. The dark shaded yellow, purple, and green regions depict 1σ regions around the retrieved cloud, gas, and Rayleigh opacity levels, respectively. The light shaded regions of the same colors depict the 2σ intervals.
Download figure:
Standard image High-resolution image

This analysis demonstrates that for observations with lower S/N (∼5), single cloud deck parameterizations are preferred when performing retrievals (even when two cloud decks are present). We also verify that all the models fail to capture any of the atmospheric properties such as molecular abundances directly with this data quality. This remains true even at higher resolutions of 140. However, inspection of the retrieved photon attenuation maps for observations with lower S/N can be informative regarding the positions of the cloud opacity and therefore the planet's climate.

5.2. Validity of Wavelength-independent Cloud Properties in Our Retrieval Models

We ignore any wavelength dependence while parameterizing the cloud optical properties for all of our retrieval models. However, when modeling our simulated data we consider full wavelength-dependent cloud optical properties from Virga. This begs the question of whether or not additional wavelength-dependent cloud complexity would be needed for our S/N and wavelength range parameter space. For the cases considered here, though, this additional complexity is not necessary because of (1) the wavelength region explored and (2) the optical properties of the cloud species explored.

The cloud optical properties of eps Eri b and 47 UMa b show negligible wavelength dependence. This is because the optical properties of H₂O and NH₃ are not strongly wavelength-dependent within the wavelength range of our focus (0.3–1 μm). Therefore, this assumption is strongly dependent on the particular species explored. For planets that are hotter than eps Eri b and 47 UMa b, condensation of sulfur-based species such as MnS, Na₂S, and ZnS may occur. These species show strong wavelength dependence in their scattering properties (e.g., Querry 1987). Although our hottest target, HD 62509 b, has a cloud deck dominated by sulfur species (Na₂S), the overall opacity is dominated by Rayleigh and molecular contributions. That is because the cloud deck of HD 62509 b is far too low in altitude (high in pressure). Therefore the necessity of wavelength-dependent properties is not warranted.

Lastly, in addition to condensates, hazes can add an additional wavelength dependence. In particular, sulfur hazes, which strongly absorb light toward 0.3 μm, can create spectra with a positive slope (Gao et al. 2017) that would require the consideration of a wavelength-dependent cloud retrieval. A Jupiter-like, wavelength-dependent haze has also been retrieved by using a simple parameterization (Lacy et al. 2018). These specific cases are beyond the scope of this analysis, but could add an additional level of complexity to the parameterizations explored here.

5.3. Additional Uncertainty Caused by Unknown Phase Angle

Throughout our analysis the phase angle of our target planets has been kept to zero. Phase, however, changes the albedo spectrum of a planet significantly because the scattering properties of the atmosphere are phase-dependent. Therefore, an unknown phase can lead to additional uncertainty in the retrieval analysis. This effect has been the subject of previous exploration. In particular, there exists a known degeneracy between the phase angle and the radius of the planet (Nayak et al. 2017). Nayak et al. (2017) showed that an unknown phase angle does not lead to a significant change in the accuracy of retrieved molecular abundances and cloud structure when compared to the case of known phase angle. However, when retrieving on contrast (relative planet-to-star ratio) as opposed to albedo, the unknown phase angle does introduce significant uncertainty in the radius retrieval of the planet compared to the case where phase angle is known (Nayak et al. 2017).

To test the sensitivity to phase angle in this analysis, we performed a retrieval on the albedo spectrum of 47 UMa b simulated at a phase angle of 90° with our Case 3 retrieval model. We assume a uniform prior for the phase angle between 60° and 120°. Similar to Nayak et al. (2017), we find that the precision and accuracy of the retrieved molecular abundances are unchanged relative to the case of zero phase, within 1σ. We do find significant bimodality in the retrieved cloud solution of (1) the pressure level of the optically thick (high pressure) cloud deck, (2) the gravity, and (3) the phase angle. Instead of retrieving a single-peaked posterior at 30 bar, a double-peaked posterior solution of 30 bar and 0.3 bar is retrieved. Therefore, the combination of unknown gravity and phase angle will impede accurate and precise determination of cloud properties. However, given that the solution is strongly bimodal, inferences could be made regarding the most likely physical scenario.

6. Conclusions and Future Work

We have performed retrievals on reflected light albedo spectra for three high-priority cool giant targets for future space-based optical high-contrast imaging and spectroscopic missions such as HabEx and LUVOIR. We have chosen planets with three different estimated effective temperatures of 135, 217, and 533 K. This wide range of effective temperatures helps to explore retrievals of various possible cloud structure scenarios for cool giants. Albedo spectra for these planets were calculated using the spectroscopy modeling code PICASO and robust cloud calculation model Virga. We used the modeled albedo spectra to simulate mock observation spectra with a constant spectral resolution of 40 and an S/N of 20. Here we briefly discuss the key aspects and results of our retrieval analysis.

1.
Requisite cloud complexity is highly sensitive to the relative position of the molecular, cloud, and Rayleigh opacities. The additional complexity of a second cloud deck, for example, is only favored (according to Bayes factor analysis) when the region with the greatest contribution to molecular opacity is between the two cloud decks. Otherwise, the atmosphere can be simply parameterized with a larger, single deck.
2.
Box model parameterizations for cloud opacity result in abundance measurements that are largely overestimated. Therefore, exponential cloud opacity parameterizations (e.g., at least Case 2) should be used instead, even at low S/N ∼ 5 observations.
3.
Although single scattering and asymmetry of the cloud deck change with altitude, a two-valued model for these scattering properties never retrieves a more accurate solution than a single-valued model (i.e., Case 3 retrieves the scattering properties more accurately than Case 4). This conclusion, however, might only apply to the planets explored here (dominated by NH₃ and H₂O clouds). Planets with two cloud decks composed of condensates or hazes with drastically different optical properties might warrant additional altitude-based complexity.
4.
Allowing for an altitude dependence in the H₂O mixing ratio profile, in order to detect H₂O depletion above the cloud deck due to water condensation for the case of 47 UMa b, slightly improves the precision and accuracy of the abundances. However, this seemingly "better" solution was weakly rejected over an identical retrieval without altitude dependence. Therefore, although the fit appears better (i.e., increased precision and accuracy with respect to the 1σ constraint interval), the additional complexity is not statistically favored for this data quality.
5.
We find that even at very low S/N = 5, low R = 40 (0.3–1 μm), inferences can be made with respect to the position of the cloud deck without attaining accurate information regarding the abundances of molecular species. In accordance with other works (Lupu et al. 2016; Nayak et al. 2017; Hu 2019) we are unable to attain precisely constrained molecular abundances with this data quality. However, we are able to retrieve a photon attenuation map of the expected opacity contribution of Rayleigh scattering, cloud scattering, and molecular absorption. This places a limit on the position of the bottom of the cloud deck. This suggests that very coarse, zeroth-order, temperature information could be attained by combining the equilibrium temperature of the planet with knowledge of condensation curves.
6.
Lastly, we show that the cloud structure and molecular mixing ratios of the planets can be accurately and precisely retrieved with a ±50% uncertainty in the gravity of the planets. However, it is not possible to improve the gravity constraint beyond this value.

We return to our initial questions posed. (1) Users' choice of cloud and atmospheric parameterization strongly affects the precision and accuracy of the resultant abundances and cloud structure. (2) The specific location of the cloud deck, with respect to the location of the optically thick molecular opacity, dictates whether or not accurate cloud structure information can be retrieved, though this information will not be known a priori. (3) Precise gravity information is very difficult to retrieve with the quality of simulated data used here, but atmospheric characterization with reflected light is possible even with large uncertainties in planet gravity. (4) Lastly, even reflected light spectroscopy with low S/N = 5 and low R = 40 from 0.3–1 μm can give insights into the position of the planet's cloud deck.

S.M. would like to thank the S. N. Bose Scholar's program by Indo-US Science and Technology Forum (IUSSTF) for funding his visit to the Department of Astronomy and Astrophysics, UC Santa Cruz through the S. N. Bose Scholarship. M.M. acknowledges the support of the Nancy Grace Roman Science Investigation Team program. The authors would like to thank the exoplanet group at UC Santa Cruz, especially Jonathan Fortney, for all the computational support and resources used in this work. The authors would also like to thank Ryan MacDonald for insightful discussions and the anonymous referee for their suggestions, which helped in improving the manuscript.

Software: PICASO (Batalha et al. 2019), Virga (Batalha et al. 2020a), DYNESTY (Speagle 2020), 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), PySynphot (STScI Development Team 2013), sqlite3, matplotlib (Hunter 2007), PyMieScatt (Sumlin et al. 2018).

Cloud Parameterizations and their Effect on Retrievals of Exoplanet Reflection Spectroscopy

Article metrics

Permissions

Share this article

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction

1.1. Previous Parameterizations Used in Reflected Light Retrievals

1.1.1. Lupu et al. (2016) and Nayak et al. (2017)

1.1.2. Feng et al. (2018)

1.1.3. Hu (2019) and Damiano & Hu (2020)

1.1.4. Irwin et al. (2008) and Barstow et al. (2014)

1.2. This Analysis

2. Modeling Reflected Spectra

2.1. Modeling the Planet Atmospheres

2.1.1. The Temperature–Pressure Profile

2.1.2. Atmospheric Chemistry

2.1.3. Clouds with Virga

2.1.4. Reflected Spectra with PICASO

2.2. The Retrieval Setup

2.2.1. The Box Cloud Model (Case 1)

2.2.2. Single Cloud Profile Model (Case 2)

2.2.3. Double Cloud Profile Model (Case 3)

2.2.4. Double Cloud Profile Model with Two-valued g0 and w0 (Case 4)

3. Dynamic Nested Sampling

4. Results

4.1. 47 UMa b

4.1.1. Comparison of Retrieval Parameterizations

4.1.2. 47 UMa b Case 3 with a Depleted Water Profile

4.1.3. Retrieval at Higher Spectral Resolution (R = 140)

4.2. eps Eri b

4.3. HD 62509 b

5. Discussion

5.1. Required Cloud Complexity

5.1.1. Ability to Constrain Gravity

5.1.2. Effect of Assumed Signal-to-noise Ratio

5.2. Validity of Wavelength-independent Cloud Properties in Our Retrieval Models

5.3. Additional Uncertainty Caused by Unknown Phase Angle

6. Conclusions and Future Work

Footnotes

2.1.3. Clouds with `Virga`

2.1.4. Reflected Spectra with `PICASO`

2.2.4. Double Cloud Profile Model with Two-valued g₀ and w₀ (Case 4)