Characterizing Earth Analogs in Reflected Light: Atmospheric Retrieval Studies for Future Space Telescopes

Our three-dimensional albedo model (see also Cahoy et al. 2010) divides a world into a number of plane-parallel facets with coordinates of longitude (ζ) and colatitude (η), with the former referenced from the subobserver location and the latter ranging from zero at the northern pole to π at the southern pole. A single facet has downwelling incident stellar radiation from a zenith angle ${\mu }_{{\rm{s}}}=\cos {\theta }_{{\rm{s}}}=\sin \eta \cos (\zeta -\alpha )$, where, as earlier, α is the phase angle. The facet reflects to the observer in a direction whose zenith angle is given by ${\mu }_{{\rm{o}}}=\cos {\theta }_{{\rm{o}}}\,=\sin \eta \cos \zeta$. Note that, at full phase (where the geometric albedo is defined), the observer and source are colinear such that ${\mu }_{{\rm{o}}}={\mu }_{{\rm{s}}}$ for all facets.

The atmosphere above each facet is divided into a set of pressure levels, and we perform a radiative transfer calculation to determine the emergent intensity. With the intensities calculated for an entire visible hemisphere, we follow the methods outlined by Horak (1950) and Horak & Little (1965) to perform integration using Chebychev–Gauss quadrature, thus producing the reflectance value at a given wavelength. We repeat this procedure at each of the wavelength points within a specified range to build up a reflectance spectrum.

Taking I(τ, μ, ϕ) to be the wavelength-dependent intensity at optical depth τ in a direction determined by the zenith and azimuth angles μ and ϕ, we ultimately need to determine the emergent intensity from each facet in the direction of the observer, I(τ = 0, μ_o, ϕ_o). Thus, for each facet, we must solve the one-dimensional, plane-parallel radiative transfer equation,

$$\begin{eqnarray}&&\mu \displaystyle \frac{{dI}}{d\tau }=I(\tau ,\mu ,\phi )-S(\tau ,\mu ,\phi ),\end{eqnarray} \tag{ 2 }$$

where S is the wavelength-dependent source function. Following Meador & Weaver (1980), Toon et al. (1989), and Marley & Robinson (2015), the source function is

$$\begin{eqnarray}&&S(\tau ,\mu ,\phi )=\displaystyle \frac{\bar{\omega }}{4\pi }{F}_{{\rm{s}}}\cdot p(\tau ,\mu ,\phi ,-{\mu }_{{\rm{s}}},{\phi }_{{\rm{s}}})\cdot {e}^{-\tau /{\mu }_{\odot }}\\ &&\quad +\,\bar{\omega }{\displaystyle \int }_{0}^{2\pi }d\phi ^{\prime} {\displaystyle \int }_{-1}^{1}\displaystyle \frac{d\mu ^{\prime} }{4\pi }\cdot I(\tau ,\mu ^{\prime} ,\phi ^{\prime} )\cdot p(\tau ,\mu ,\phi ,\mu ^{\prime} ,\phi ^{\prime} ),\end{eqnarray} \tag{ 3 }$$

where $\bar{\omega }$ is the single-scattering albedo, F_s is the incoming stellar flux at the top of the atmosphere (which we normalize to unity so that emergent intensities correspond to reflectivities), ϕ_s is the stellar azimuth angle, and p is the scattering phase function. Note that our source function does not include an emission term, since we are not computing thermal spectra. Recall that the first term in Equation (3) describes directly scattered radiation from the direct solar beam, while the final term describes diffusely scattered radiation from the ($\mu ^{\prime} ,\phi ^{\prime}$) direction scattering into the (μ, ϕ) direction.

Like most standard tools for solving the radiative transfer equation, we separate treatments of directly scattered radiation from diffusely scattered radiation, and, for both, it is convenient to express the scattering phase function in terms of a unique scattering angle, Θ. As single-scattered radiation typically has more distinct forward- and backward-scattering features, we choose to represent the scattering phase function for the direct beam with a two-term Henyey–Greenstein (TTHG) phase function (Kattawar 1975),

$$\begin{eqnarray}&&{p}_{\mathrm{TTHG}}({\rm{\Theta }})={{fp}}_{\mathrm{HG}}({g}_{{\rm{f}}},{\rm{\Theta }})+(1-f){p}_{\mathrm{HG}}({g}_{{\rm{b}}},{\rm{\Theta }})\,,\end{eqnarray} \tag{ 4 }$$

where p_HG is the Henyey–Greenstein (HG) phase function with

$$\begin{eqnarray}&&{p}_{\mathrm{HG}}=\displaystyle \frac{1}{4\pi }\displaystyle \frac{1-{\bar{g}}^{2}}{{(1+{\bar{g}}^{2}-2\bar{g}\cos {\rm{\Theta }})}^{3/2}}.\end{eqnarray} \tag{ 5 }$$

Recall that a TTHG phase function can represent both forward- and backward-scattering peaks, while the (one-term) HG phase function only has one peak (typically in the forward direction). In the previous expressions, $\bar{g}$ is the asymmetry parameter, f is the forward/backward scatter fraction, g_f is the asymmetry parameter for the forward-scattered portion of the TTHG, and g_b is the asymmetry parameter for the backward-scattered portion of the TTHG. For the forward- and backward-scattering portions of the TTHG phase function, we use ${g}_{{\rm{f}}}=\bar{g}$, ${g}_{{\rm{b}}}=-\bar{g}/2$, and $f=1-{g}_{{\rm{b}}}^{2}$. Substituting these into the TTHG phase function expression yields

$$\begin{eqnarray}&&{p}_{\mathrm{TTHG}}=\displaystyle \frac{{\bar{g}}^{2}}{4}{p}_{\mathrm{HG}}(-\bar{g}/2,{\rm{\Theta }})+\left(1-\displaystyle \frac{{\bar{g}}^{2}}{4}\right){p}_{\mathrm{HG}}(\bar{g},{\rm{\Theta }}).\end{eqnarray} \tag{ 6 }$$

For radiation that is single-scattered from the solar beam to the observer, the scattering geometry is fixed by the planetary phase angle such that Θ = π − α. Our choice of parameters, and their relation to $\bar{g}$, in the TTHG was designed by Cahoy et al. (2010) to roughly reproduce the phase function of liquid-water clouds. This parameterization, however, is different from that proposed by Kattawar (1975). We do not expect our results to be sensitive to the details of a particular phase function treatment, as Lupu et al. (2016) showed that scattered-light retrievals struggle to constrain phase function parameters.

We adopt a standard two-stream approach to solving the radiative transfer equation (Meador & Weaver 1980). In this case, the diffusely scattered component of the source function is azimuthally averaged. Combined with our representation of the directly scattered component, we have

$$\begin{eqnarray}S(\tau ,\mu ,{\mu }_{{\rm{s}}},\alpha ) & = & \displaystyle \frac{\bar{\omega }}{4\pi }{F}_{{\rm{s}}}\cdot {p}_{\mathrm{TTHG}}(\mu ,-{\mu }_{{\rm{s}}})\cdot {e}^{-\tau /{\mu }_{{\rm{s}}}}\\ & & +\displaystyle \frac{\bar{\omega }}{4\pi }{\displaystyle \int }_{-1}^{1}I(\tau ,\mu ^{\prime} )p(\mu ,\mu ^{\prime} )d\mu ^{\prime} ,\end{eqnarray} \tag{ 7 }$$

where the azimuth-averaged phase functions are given by

$$\begin{eqnarray}&&p(\mu ,\mu ^{\prime} )=\displaystyle \frac{1}{2\pi }{\int }_{0}^{2\pi }p(\mu ,\phi ,\mu ^{\prime} ,\phi )d\phi .\end{eqnarray} \tag{ 8 }$$

We represent the azimuth-averaged scattering phase functions as a series of Legendre polynomials, P_l(μ), expanded to order M with

$$\begin{eqnarray}&&p(\mu ,\mu ^{\prime} )=\displaystyle \sum _{l=0}^{M}{g}_{l}{P}_{l}(\mu ){P}_{l}(\mu ^{\prime} ),\end{eqnarray} \tag{ 9 }$$

where the phase function moments, g_l, are defined according to

$$\begin{eqnarray}&&{g}_{l}=\displaystyle \frac{2l+1}{2}{\int }_{-1}^{1}p(\cos {\rm{\Theta }}){P}_{l}(\cos {\rm{\Theta }})d\cos {\rm{\Theta }}.\end{eqnarray} \tag{ 10 }$$

The first moment of the phase function is related to the asymmetry parameter, with $\bar{g}={g}_{1}/3$. We use a second-order expansion of the phase function, giving

$$\begin{eqnarray}&&p(\mu ,\mu ^{\prime} )=1+3\bar{g}\mu \mu ^{\prime} +\displaystyle \frac{{g}_{2}}{2}(3{\mu }^{2}-1)(3\mu {{\prime} }^{2}-1).\end{eqnarray} \tag{ 11 }$$

In a given atmospheric layer of our albedo model, the optical depth is the sum of the scattering optical depth and the absorption optical depth, τ = τ_scat + τ_abs. The scattering optical depth has contributions from Rayleigh scattering and clouds, so that ${\tau }_{\mathrm{scat}}={\tau }_{\mathrm{Ray}}+{\bar{\omega }}_{\mathrm{cld}}{\tau }_{\mathrm{cld}}$, where ${\bar{\omega }}_{\mathrm{cld}}$ is the cloud single-scattering albedo. The single-scattering albedo for a layer is then $\bar{\omega }={\tau }_{\mathrm{scat}}/\tau$. We determine the asymmetry parameter, $\bar{g}$, with an optical depth weighting on the Rayleigh-scattering asymmetry parameter (which is zero) and the cloud-scattering asymmetry parameter, yielding $\bar{g}\,={\bar{g}}_{\mathrm{cld}}({\tau }_{\mathrm{cld}}/{\tau }_{\mathrm{scat}})$. When representing the second moment of the phase function, we use ${g}_{2}=\tfrac{1}{2}({\tau }_{\mathrm{Ray}}/{\tau }_{\mathrm{scat}})$ so that g₂ tends toward the appropriate value for Rayleigh scattering (i.e., 1/2; Hansen & Travis 1974) when the Rayleigh-scattering optical depth dominates the scattering optical depth.

As compared to prior investigations that have used the Cahoy et al. (2010) albedo tool (e.g., Lupu et al. 2016; Nayak et al. 2017), we have updated the model to include an optional isotropically reflecting (Lambertian) lower boundary (mimicking a planetary surface) and added pressure-dependent absorption due to ${{\rm{H}}}_{2}{\rm{O}}$, O₃, O₂, and CO₂. As in previous studies, CH₄ remains a radiatively active species in the model. We also include Rayleigh scattering from ${{\rm{H}}}_{2}{\rm{O}}$, O₂, CO₂, and N₂ (in addition to H₂ and He from previous studies). As in Lupu et al. (2016), we allow for an extended gray-scattering cloud in our atmospheres.

In Lupu et al. (2016), their two-layer cloud model atmosphere includes a deeper, optically thick cloud deck that essentially acts as a reflective surface. Unlike the gas giants within that study, though, terrestrial planets have a solid surface we can probe. We characterize our isotropically reflecting lower boundary using a spherical albedo for the planetary surface, ${A}_{{\rm{s}}}$, which represents the specific power in scattered, outgoing radiation compared to that in incident radiation. For this study, we simply adopt gray surface albedo values, which reduces complexity and computation time. For the inhomogeneous surface of a realistic Earth, featuring oceans and continents, the surface albedo is wavelength-dependent, and we hope to investigate the significance of such surfaces in future work.

We undertook a test to check our reflective lower-boundary condition in the limit of a transparent atmosphere. Without atmospheric absorption or scattering, our assumption of a Lambertian surface would imply that the reflectivity (or phase function) determined by our albedo code should follow the analytic Lambert phase function,

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{{\rm{L}}}(\alpha )=\displaystyle \frac{\sin \alpha +(\pi -\alpha )\cos \alpha }{\pi }.\end{eqnarray} \tag{ 12 }$$

Figure 1 compares the model phase function with the analytic phase function and shows complete agreement, confirming that our treatment of the surface is correct.

Zoom In Zoom Out Reset image size

Previous work featuring the albedo model adopted here used a predefined atmospheric pressure grid. To accommodate the finite surface pressures of rocky planets, as well as the various combinations of cloud parameters our retrievals will explore, we instead establish an adaptive method of determining the pressure grid. Here we divide the atmosphere into a pressure grid of N_level, bounded by P = P_top at the top of the atmosphere and P = P₀ at the surface. In a cloud-free scenario, we simply divide the atmosphere evenly in log-P space.

For our simulations that include a single cloud deck, we adaptively determine the pressure value at each level depending on the location, thickness, and optical depth of the cloud. The quantities that define the cloud deck are p_t, the cloud-top pressure; dp, the atmospheric pressure across the cloud; and τ, the cloud optical depth. We begin by assigning a number of layers to the cloud, imposing two conditions: (1) there should be at least three model pressure layers to each atmospheric pressure scale height (perH = 3), and (2) the cloud optical depth in a layer must remain below at most 5 (maxtau = 5). This allows us to avoid any one layer spanning a large extent within the atmosphere and also avoids cloud layers that have extremely large scattering optical depths.

When beginning our gridding process, we propose an initial number of cloud layers, ${N}_{{\rm{c}}}={\mathtt{perH}}\times {\mathtt{numH}}$, where ${\mathtt{numH}}=\mathrm{ln}\tfrac{{p}_{{\rm{t}}}+{dp}}{{p}_{{\rm{t}}}}$ is the number of e-folding distances through the cloud (serving as a proxy for scale height). The aerosol optical depth for each pressure layer within the cloud would then simply be ${\rm{\Delta }}\tau =\tfrac{\tau }{{N}_{{\rm{c}}}}$. However, if ${\rm{\Delta }}\tau \gt {\mathtt{maxtau}}$, we adjust the cloud resolution by increasing N_c by a factor of $\tfrac{{\rm{\Delta }}\tau }{{\mathtt{maxtau}}}$ and then round up to the nearest integer. In other words, we increase the resolution of the pressure grid through the cloud until the layer optical depth is under maxtau. We determine successive pressure level values through the cloud with $p[i]=p[i-1]+{\rm{\Delta }}\mathrm{ln}p$, where ${\rm{\Delta }}\mathrm{ln}p=\tfrac{\mathrm{ln}({p}_{{\rm{t}}}+{dp})-\mathrm{ln}{p}_{{\rm{t}}}}{{N}_{{\rm{c}}}}$, starting from the top of the cloud. We divide the remaining ${N}_{\mathrm{level}}-{N}_{{\rm{c}}}$ levels in uniform $\mathrm{ln}p$ space on either side of the cloud, weighted by the number of pressure scale heights above (${N}_{{\rm{t}}}$) and below (${N}_{{\rm{b}}}$) the cloud. Figure 2 visualizes the three portions of the atmosphere.

Zoom In Zoom Out Reset image size

For simplicity, we assume an isothermal atmosphere (at T = 250 K), as temperature has little effect on the reflected-light spectrum (Robinson 2017). Pressure, however, has a strong impact on molecular opacities, as seen in Figure 3. We incorporated high-resolution pressure-dependent opacities for all molecules in our atmosphere. The absorption opacities are generated line-by-line from the HITRAN2012 line list (Rothman et al. 2013) for seven orders of magnitude in pressure (10⁻⁵–10² bar) at T, spanning our entire wavelength range at <1 cm⁻¹ resolution. Figure 3 also illustrates how absorption features of ${{\rm{H}}}_{2}{\rm{O}}$, O₂, and O₃ change when in an atmosphere of 1 bar versus one of 10 bars.

Zoom In Zoom Out Reset image size

We interpolate our high-resolution opacity tables to the slightly lower resolution of the forward model in order to maintain short model runtimes while not affecting the accuracy of the output spectra. For each model layer, we interpolate over the opacities from our table given the pressure. The chemical abundances in our forward-model atmosphere are constant as a function of pressure, and we also adopt a uniform acceleration due to gravity.

We have also added an option to include partial cloudiness across a planetary disk whose fractional coverage is described by f_c. To mimic partial cloudiness as we see on Earth, we call the forward model twice. We use the same set of atmospheric and planetary parameters for both calls, except for the cloud optical depth. "Cloudy" is the call that has a nonzero cloud optical depth, while "cloud-free" is the call where we set the cloud optical depth to zero. Each call returns a geometric albedo spectrum, and we combine the two sets with the fractional cloudiness parameter such that the combined spectrum follows ${f}_{c}\times \mathrm{cloudy}+(1-{f}_{c})\times \mathrm{cloud}\mbox{--}\mathrm{free}$.

The generalized three-dimensional albedo model described above can simulate reflected-light spectra of a large diversity of planet types, spanning solid-surfaced worlds to gas giants with a variety of prescribed atmospheric compositions. For the present study, however, we choose to focus on Earth-like worlds, which are described in detail below. Thus, we define a set of fiducial model input parameters that are designed to mimic Earth and thereby enable us to generate simulated observational data sets for an Earth twin.

Table 1 summarizes the fiducial model parameter values adopted for our Earth twin. Also shown are the priors for these parameters, which we use when performing retrieval analyses. For an Earth-like setup, the surface atmospheric pressure is P₀ = 1 bar, and we adopt a surface albedo of A_s = 0.05, which is representative of mostly ocean-covered surface. We adopt a uniform acceleration due to gravity of g = 9.8 m s⁻² and set the planetary radius to R_⊕. For convenience, we sometime refer to these four variables (P₀, A_s, g, and R_p) as the bulk planetary and atmospheric parameters.

Table 1.  List of the 11 Retrieved Parameters in the Complete Cloudy Model, Their Descriptions, Fiducial Input Values, and Corresponding Priors

Parameter	Description	Input	Prior
$\mathrm{log}{P}_{0}$ (bar)	Surface pressure	$\mathrm{log}(1)$	[−2,2]
${\mathrm{logH}}_{2}{\rm{O}}$	Water vapor mixing ratio	$\mathrm{log}(3\times {10}^{-3})$	[−8,−1]
${\mathrm{logO}}_{3}$	Ozone mixing ratio	$\mathrm{log}(7\times {10}^{-7})$	[−10,−1]
${\mathrm{logO}}_{2}$	Molecular oxygen mixing ratio	$\mathrm{log}(0.21)$	[−10,0]
Rp (${R}_{\oplus }$)	Planet radius	1	[0.5, 12]
$\mathrm{log}g$ (m s−2)	Surface gravity	$\mathrm{log}(9.8)$	[0,2]
$\mathrm{log}{A}_{s}$	Surface albedo	$\mathrm{log}(0.05)$	[−2, 0]
$\mathrm{log}{p}_{t}$ (bar)	Cloud-top pressure	$\mathrm{log}(0.6)$	[−2,2]
$\mathrm{log}{dp}$ (bar)	Cloud thickness	$\mathrm{log}(0.1)$	[−3,2]
$\mathrm{log}\tau$	Cloud optical depth	$\mathrm{log}(10)$	[−2,2]
$\mathrm{log}{f}_{c}$	Cloudiness fraction	$\mathrm{log}(0.5)$	[−3,0]

Download table as:  ASCIITypeset image

We focus on molecular absorption due to ${{\rm{H}}}_{2}{\rm{O}}$, O₃, and O₂. While our albedo model includes opacities from CH₄ and CO₂ as well, we omit these two species, as the reflected-light spectrum of Earth in the visible contains no strong features for these molecules. The input values for the molecular abundances (or volume mixing ratios) are ${{\rm{H}}}_{2}{\rm{O}}$ = 3 × 10⁻³, O₃ = 7 × 10⁻⁷, and O₂ = 0.21. These abundance values are based on column-weighted averages from a standard Earth model atmosphere with vertically varying gas mixing ratios (McClatchey et al. 1972). The primary Rayleigh scatterer and background gas in our fiducial model is N₂, whose abundance makes up the remainder of the atmosphere after all other gases are accounted for (i.e., roughly 0.79). Rayleigh scattering is treated according to Hansen & Travis (1974) with constants to describe the scattering properties of N₂, O₂, and ${{\rm{H}}}_{2}{\rm{O}}$ from Allen & Cox (2000). We do not include polarization or Raman-scattering effects.

Our cloud model was designed to be minimally parametric while still enabling us to sufficiently reproduce realistic spectra of Earth. Our single-layer gray ${{\rm{H}}}_{2}{\rm{O}}$ cloud has $\bar{\omega }=1$ and $\bar{g}=0.85$, which are characteristic of water clouds across the visible range. These two parameters were fixed to minimize retrieval model complexity, as we believe that water is the most likely condensate for worlds in the habitable zone. Nevertheless, future studies may not wish to assume values of $\bar{\omega }$ and $\bar{g}$ a priori. Cloud-top pressure (p_t) and fractional coverage (f_c) are set at 0.6 bar and 50%, respectively, which are roughly consistent with observations of optically thick cloud coverage on Earth (Stubenrauch et al. 2013). Cloud thickness (dp) and optical depth (τ) were set to 0.1 bar and 10, respectively, based on results from the MODIS instrument (http://modis-atmos.gsfc.nasa.gov) used in Robinson et al. (2011).

With fiducial values chosen, we validate our forward model against a simulated high-resolution disk-integrated spectrum of Earth at full phase, as shown in Figure 4. The comparison spectrum is produced by the NASA Astrobiology Institute's Virtual Planetary Laboratory (VPL) sophisticated 3D, line-by-line, multiple-scattering spectral Earth model (Robinson et al. 2011). The Robinson et al. (2011) tool can simulate images and disk-integrated spectra of Earth from the ultraviolet to the infrared. It has been validated against observations at visible wavelengths taken by NASA's EPOXI (Robinson et al. 2011) and LCROSS missions (Robinson et al. 2014).

Zoom In Zoom Out Reset image size

Features of the Robinson et al. (2011) model include Rayleigh scattering due to air molecules, realistic patchy clouds, and gas absorption from a variety of molecules, including ${{\rm{H}}}_{2}{\rm{O}}$, CO₂, O₂, O₃, and CH₄. Surface coverage of different land types (e.g., forest, desert) is informed by satellite data, and water surfaces incorporate specular reflectance of sunlight. A grid of thousands of surface pixels is nested beneath a grid of 48 independent atmospheric pixels, all of equal area. For each surface pixel, properties from the overlying atmospheric pixels are used as inputs to a full-physics, plane-parallel radiative transfer solver: the Spectral Mapping Atmospheric Radiative Transfer (SMART) model (Meadows & Crisp 1996). Intensities from this solver are integrated over the pixels with respect to solid angle, thereby returning a disk-integrated spectrum.

The sophistication of the Robinson et al. (2011) model makes it unsuitable to retrieval studies, however, as model runtimes are measured in weeks for the highest-complexity scenarios. This, in part, justifies our adoption of a minimally parametric albedo model (whose runtime is measured in seconds). Furthermore, as in Figure 4, our efficient albedo model reproduces all of the key features of the Robinson et al. (2011) model. The most notable differences—that the efficient model, as compared to the Robinson et al. (2011) model, is more reflective in the blue and less reflective in the red—are simply due to our adoption of a gray surface albedo. Land and plants, which cover roughly 29% of Earth's surface, are generally more reflective in the red than in the blue. Figure 4 also compares a spectrum from our forward model against a spectrum of a partially clouded ocean planet generated with the Robinson et al. (2011) model. This ocean world is identical to Earth except for the fact that its surface is covered entirely by an ocean, with no land present. The surface albedo in the ocean model is gray beyond 500 nm; shortward of this, the reflectivity increases, likely leading to the discrepancy in our comparison at the bluest wavelengths. Still, with a more accurate match to a planet that has a nearly gray albedo through the visible, we consider our assumption of gray surface albedo to be the main reason for the discrepancies when compared to the Robinson et al. (2011) realistic model.

Finally, in our albedo model, we set 100 facets for the visible hemisphere and calculate a high-resolution geometric albedo spectrum at 1000 wavelength points between 0.35 and 1.05 μm. Like Lupu et al. (2016), we only consider a planet at full phase (α = 0°). While direct-imaging missions will not obtain observations of exoplanets at full phase, this assumption makes little difference for our results, as we are not computing integration times and only work in S/N space. Additionally, as Nayak et al. (2017) followed up Lupu et al. (2016) by retrieving phase information from giant planets in reflected light, we anticipate performing a similar expansion in the future. Our forward model has 61 pressure levels in an isothermal atmosphere of 250 K, bounded below by a reflective surface. The top of the atmosphere is set at P_top = 10⁻⁴ bars.

We convert a high-resolution geometric albedo spectrum to a synthetic planet-to-star flux ratio spectrum given the resolution of an instrument and a noise model. We then apply a Bayesian inference tool on the synthetic data set to sample the posterior probability distributions of the forward-model input parameters. To perform Bayesian parameter estimation, we utilize the open-source affine-invariant MCMC ensemble sampler emcee (Goodman & Weare 2010; Foreman-Mackey et al. 2013). Ensemble refers to the use of many chains, or walkers, to traverse parameter space; as a massively parallelized algorithm, it is computationally efficient. Affine invariance refers to the invariant performance under linear transformations of parameter space, enabling the algorithm to be insensitive to parameter covariances (Foreman-Mackey et al. 2013). With a cloudy retrieval, we can expect complex correlations that a sampler should be able to reveal. As it is more agnostic to the shape of the posterior, we choose emcee following Nayak et al. (2017) over Multinest, another sampler Lupu et al. (2016) considered that yielded consistent results. The albedo model is coded in Fortran; we convert it into a Python-callable library with the F2PY package. Each call to the forward model takes approximately 10 s of clock time on an eight-core processor. To visualize the MCMC results, we utilize the corner-plotting package developed by Foreman-Mackey (2016).

Table 1 lists the priors for our parameters. We offer a generous range on the molecular abundances; we allow O₂ in particular to dominate the atmospheric composition. Our choice of radius range (0.5–12 ${R}_{\oplus }$) reflects the range of planetary sizes from Mars to Jupiter. Also, when performing retrievals, we impose two limiting conditions to maintain physical scenarios. First, we limit the mixing ratio of N₂, ${f}_{{{\rm{N}}}_{2}}\,=1-\sum \,(\mathrm{gas}\ \mathrm{abundances})$, to be between 0 and 1. Second, for the cloud pressure terms, we reject any drawn value that does not satisfy ${10}^{{p}_{{\rm{t}}}}+{10}^{{dp}}\lt {10}^{{P}_{0}}$ (i.e., that the cloud base cannot extend below the bottom of the atmosphere). Note that for the purposes of the retrieval, we consider pressures in log space.

We simulate noise in our observations following the expressions given in Robinson et al. (2016). For simplicity, we include only read noise and dark current, as Robinson et al. (2016) showed that detector noise will be the dominant noise source in WFIRST-type spectral observations of exoplanets. Detector and instrument parameters for the HabEx and LUVOIR concepts are only loosely defined, and advances in detector technologies for these missions may move observations out of the detector-noise-dominated regime. In the detector-noise-dominated regime, the S/N is simply

$$\begin{eqnarray}&&{\rm{S}}/{\rm{N}}=\displaystyle \frac{{c}_{p}\times {t}_{\mathrm{int}}}{\sqrt{({c}_{d}+{c}_{r})\times {t}_{\mathrm{int}}}}\,\end{eqnarray} \tag{ 13 }$$

,where t_int is the integration time, c_p is the planet count rate, c_d is the dark-noise count rate, and c_r is the read-noise count rate. More rigorously, it can be shown that, at constant spectral resolution, ${\rm{S}}/{\rm{N}}\propto q{ \mathcal T }{A}_{g}{\rm{\Phi }}(\alpha ){B}_{\lambda }\lambda$, where q is the wavelength-dependent detector quantum efficiency, ${ \mathcal T }$ is throughput, and B_λ is the host stellar specific intensity (taken here as a Planck function at the stellar effective temperature). We use a stellar temperature of 5780 K for the blackbody. When the S/N at one wavelength is specified, this scaling implies that the calculation of the S/N at other wavelengths is independent of the imaging raw contrast of the instrument. We can expect the noise at the redder end of our range to be large, as the detector quantum efficiency (taken to be appropriate for the WFIRST/CGI) rapidly decreases. Since we treat only S/N rather than modeling exposure times, the exact mix of noise sources is not relevant (so that, e.g., dark current and read noise are indistinguishable). The key relevant properties of the noise model are that it is uncorrelated between spectral channels and its magnitude only depends on wavelength via a dependence on point-spread function area (Robinson et al. 2016, their Equation (26)), which will be true for detector-limited cases but may not be true for large-aperture instruments limited by speckle noise.

For our study, we will consider multiple wavelength resolutions, R, and S/Ns. Working in S/Ns (instead of integration times) makes our investigations independent of telescope diameter, target distance, and other system-specific or observing parameters. Because the S/N is dependent on wavelength, we reference our values to be at the V band (550 nm) for all resolutions for HabEx/LUVOIR. Since the WFIRST/CGI spectrograph is currently planned to only extend to 600 nm at the blue end, we opt to reference our WFIRST S/Ns to this wavelength. Unlike previous studies (Lupu et al. 2016; Nayak et al. 2017), our simulated WFIRST rendezvous data include two photometric points in the blue, which is consistent with current CGI designs. We set the S/N in the WFIRST filters to be equal to that at 600 nm.

Our simulation grid setup is shown in Table 2, where the spectral resolutions and S/Ns assumed for different observing scenarios are indicated. Figure 5 demonstrates the WFIRST rendezvous scenario data along with R = 70 and 140 data points (for HabEx/LUVOIR) plotted over the forward-model spectrum before noise is added. The scaling of the S/N with wavelength for the WFIRST rendezvous (normalized to unity at 600 nm), as well as our R = 70 and 140 cases (normalized to unity at 550 nm), is shown in Figure 6. The impact of the host stellar SED sets the overall shape of the S/N scaling, with additional influence from atmospheric absorption band features, as well as detector quantum efficiency effects (that have strong impacts at red wavelengths). Thus, Figure 6 can be used to translate our stated S/N to the S/N at any other wavelength (e.g., an S/N = 10 simulation has an S/N in the continuum shortward of the 950 nm water vapor band of roughly 0.3 × 10 = 3).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

When generating simulated data with a noise model, there are several options for handling the placement/sampling of the mock observational data points. Previous studies (Lupu et al. 2016; Nayak et al. 2017) have generated a single randomized data set for a given S/N. The placement of a single spectral data point is determined by randomly sampling a Gaussian distribution whose width is determined by the wavelength-dependent S/N. While this treatment can accurately simulate a single observational instance, it also runs the risk (especially at lower spectral resolution and S/N) of biasing retrieval results, as the random placement of only a small handful of spectral data points can significantly impact the outcome. Given this, it is ideal to retrieve on a large number (≳10) of simulated data sets at a given spectral resolution and S/N, where a comprehensive view of all the posteriors from the collection of instances will indicate expected telescope/instrument performance. Unfortunately, given the large number of R–S/N pairs in our study (10) and the long runtime of an individual retrieval (of order 1 week on a cluster), running ∼10 noise instances for each of our R–S/N pairs is computationally unfeasible (requiring ∼100 weeks of cluster time). Thus, we opt for an intermediate approach that maintains computational feasibility and avoids potential biases from individual noise instances. Here we run only a single noise instance at a given R–S/N pair, but we do not randomize the placement of the individual spectral points. In other words, the individual simulated spectral points are placed on the "true" planet-to-star flux ratio point and assigned error bars according to the S/N and noise model. While this approach prevents having a small handful of randomized data points biasing the retrieval results, it does lead to likely optimistic results, especially at modest S/N (i.e., S/N ∼ 10), since data point randomization is, in effect, an additional "noise" source that we are omitting. This means that the posterior distributions will usually be centered on the true values in an unrealistic fashion. However, the width and shape of the posterior covariances will be representative of real observations, so the fidelity of retrievals can be assessed. We keep this optimism in mind when discussing results in later sections; in particular, we compare the performance of retrievals on multiple noise instances of a subset of the cases we consider to the nonrandomized case in Section 5.2.

For both R = 70 and R = 140, we simulated data sets at S/N = 5, 10, 15, and 20. Table 4 lists the median and 1σ values of all retrieved parameters for each S/N at R = 70. Figure 11 shows the marginalized posterior distributions for the model parameters for all S/N cases for R = 70, plotted with the fiducial or "truth" values. Table 5 lists the median and 1σ values of all retrieved parameters for each S/N at R = 140. Figure 12 shows the posterior distributions for R = 140 for the model parameters for all S/N cases compared against their input values. Figure 14 shows the corresponding spread in fits and the median fit to the data for each S/N for both resolutions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4.  R = 70 Retrieval Results, with Median Value and 1σ Uncertainties of the Parameters

Parameter	Input	S/N = 5	S/N = 10	S/N = 15	S/N = 20
$\mathrm{log}\,{{\rm{H}}}_{2}{\rm{O}}$	−2.52	$-{5.07}_{-1.92}^{+2.34}$	$-{3.85}_{-2.60}^{+1.77}$	$-{3.12}_{-1.71}^{+0.97}$	$-{2.76}_{-0.88}^{+0.62}$
$\mathrm{log}\,{{\rm{O}}}_{3}$	−6.15	$-{7.55}_{-1.46}^{+1.49}$	$-{6.79}_{-1.81}^{+0.93}$	$-{6.37}_{-0.84}^{+0.55}$	$-{6.24}_{-0.60}^{+0.47}$
$\mathrm{log}\,{{\rm{O}}}_{2}$	−0.68	$-{5.12}_{-3.23}^{+3.25}$	$-{4.51}_{-3.61}^{+3.24}$	$-{1.86}_{-3.99}^{+1.29}$	$-{1.00}_{-1.01}^{+0.66}$
$\mathrm{log}\,{{\rm{P}}}_{0}$	0.0	${0.02}_{-0.84}^{+1.35}$	$-{0.03}_{-0.70}^{+0.87}$	${0.28}_{-0.56}^{+0.85}$	${0.25}_{-0.49}^{+0.56}$
${R}_{{\rm{p}}}$	1.0	${1.23}_{-0.58}^{+1.54}$	${1.33}_{-0.52}^{+1.23}$	${0.97}_{-0.27}^{+0.68}$	${0.98}_{-0.25}^{+0.44}$
$\mathrm{log}\,g$	0.99	${1.33}_{-0.77}^{+0.48}$	${1.48}_{-0.68}^{+0.38}$	${1.28}_{-0.66}^{+0.51}$	${1.24}_{-0.69}^{+0.55}$
$\mathrm{log}\,{A}_{s}$	−1.3	$-{0.96}_{-0.74}^{+0.58}$	$-{1.05}_{-0.59}^{+0.55}$	$-{0.70}_{-0.62}^{+0.37}$	$-{0.63}_{-0.46}^{+0.29}$
$\mathrm{log}\,{p}_{t}$	−0.22	$-{1.14}_{-0.61}^{+0.97}$	$-{1.19}_{-0.56}^{+0.93}$	$-{0.92}_{-0.71}^{+0.86}$	$-{0.94}_{-0.73}^{+0.84}$
$\mathrm{log}\,{dp}$	−1.0	$-{1.67}_{-0.92}^{+1.24}$	$-{1.71}_{-0.91}^{+1.18}$	$-{1.35}_{-1.14}^{+1.17}$	$-{1.43}_{-1.06}^{+1.11}$
$\mathrm{log}\,\tau$	1.0	${0.10}_{-1.43}^{+1.30}$	${0.21}_{-1.48}^{+1.23}$	${0.49}_{-1.66}^{+1.03}$	${0.61}_{-1.66}^{+0.93}$
$\mathrm{log}\,{f}_{c}$	−0.3	$-{1.43}_{-1.07}^{+0.99}$	$-{1.33}_{-1.12}^{+0.94}$	$-{0.93}_{-1.32}^{+0.71}$	$-{1.05}_{-1.27}^{+0.80}$

Download table as:  ASCIITypeset image

Table 5.  R = 140 Retrieval Results, with Median Value and 1σ Uncertainties of the Parameters

Parameter	Input	S/N = 5	S/N = 10	S/N = 15	S/N = 20
$\mathrm{log}\,{{\rm{H}}}_{2}{\rm{O}}$	−2.52	$-{4.56}_{-2.35}^{+2.14}$	$-{2.74}_{-1.07}^{+0.69}$	$-{2.61}_{-0.65}^{+0.47}$	$-{2.43}_{-0.56}^{+0.39}$
$\mathrm{log}\,{{\rm{O}}}_{3}$	−6.15	$-{7.36}_{-1.65}^{+1.26}$	$-{6.26}_{-0.68}^{+0.53}$	$-{6.18}_{-0.48}^{+0.42}$	$-{6.03}_{-0.48}^{+0.34}$
$\mathrm{log}\,{{\rm{O}}}_{2}$	−0.68	$-{4.45}_{-3.69}^{+3.08}$	$-{1.06}_{-1.43}^{+0.76}$	$-{0.76}_{-0.79}^{+0.51}$	$-{0.60}_{-0.59}^{+0.43}$
$\mathrm{log}\,{{\rm{P}}}_{0}$	0.0	${0.07}_{-0.84}^{+1.01}$	${0.20}_{-0.49}^{+0.72}$	${0.12}_{-0.36}^{+0.49}$	${0.07}_{-0.31}^{+0.39}$
${R}_{{\rm{p}}}$	1.0	${1.25}_{-0.52}^{+1.16}$	${1.01}_{-0.28}^{+0.60}$	${0.99}_{-0.23}^{+0.42}$	${1.05}_{-0.27}^{+0.42}$
$\mathrm{log}\,g$	0.99	${1.36}_{-0.74}^{+0.46}$	${1.31}_{-0.77}^{+0.49}$	${1.14}_{-0.65}^{+0.56}$	${1.20}_{-0.64}^{+0.50}$
$\mathrm{log}\,{A}_{s}$	−1.3	$-{0.98}_{-0.60}^{+0.54}$	$-{0.67}_{-0.50}^{+0.32}$	$-{0.68}_{-0.44}^{+0.29}$	$-{0.79}_{-0.69}^{+0.34}$
$\mathrm{log}\,{p}_{t}$	−0.22	$-{1.23}_{-0.55}^{+1.03}$	$-{0.96}_{-0.71}^{+0.80}$	$-{0.79}_{-0.82}^{+0.70}$	$-{0.66}_{-0.85}^{+0.53}$
$\mathrm{log}\,{dp}$	−1.0	$-{1.72}_{-0.91}^{+1.25}$	$-{1.43}_{-1.09}^{+1.13}$	$-{1.55}_{-1.00}^{+1.10}$	$-{1.49}_{-0.98}^{+1.00}$
$\mathrm{log}\,\tau$	1.0	${0.18}_{-1.49}^{+1.31}$	${0.50}_{-1.66}^{+1.09}$	${0.61}_{-1.61}^{+0.98}$	${0.79}_{-1.40}^{+0.87}$
$\mathrm{log}\,{f}_{c}$	−0.3	$-{1.30}_{-1.13}^{+0.93}$	$-{1.31}_{-1.21}^{+0.94}$	$-{0.99}_{-1.27}^{+0.76}$	$-{0.76}_{-1.26}^{+0.54}$

Download table as:  ASCIITypeset image

For the WFIRST rendezvous scenario, we utilized the Design Cycle 7 instrument parameters to set the locations of the two photometric points and the range and resolution of the spectrometer (R = 50; see Table 2). Because of this particular setup, we reference the S/Ns in our grid (5, 10, 15, 20) at 600 nm and assign the photometric points the same S/N as at 600 nm. Table 6 lists the median and 1σ values of all retrieved parameters for each S/N variant. Figure 13 presents the posterior distributions for the four WFIRST rendezvous variants with respect to the input values. Figure 14 shows the spread in fits and median fit to the data for each variant.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 6.  WFIRST Rendezvous Retrieval Results, with Median Value and 1σ Uncertainties of the Parameters

Parameter	Input	S/N = 5	S/N = 10	S/N = 15	S/N = 20
$\mathrm{log}\,{{\rm{H}}}_{2}{\rm{O}}$	−2.52	$-{4.94}_{-2.05}^{+2.35}$	$-{4.89}_{-2.11}^{+2.48}$	$-{4.03}_{-2.52}^{+1.87}$	$-{3.11}_{-1.71}^{+1.17}$
$\mathrm{log}\,{{\rm{P}}}_{0}$	0.0	$-{0.16}_{-0.80}^{+1.32}$	$-{0.19}_{-0.71}^{+1.03}$	${0.03}_{-0.74}^{+1.16}$	${0.45}_{-0.85}^{+1.01}$
$\mathrm{log}\,{{\rm{O}}}_{3}$	−6.15	$-{7.66}_{-1.59}^{+1.65}$	$-{7.53}_{-1.66}^{+1.54}$	$-{7.16}_{-1.66}^{+1.19}$	$-{6.80}_{-1.30}^{+0.94}$
$\mathrm{log}\,{{\rm{O}}}_{2}$	−0.68	$-{5.05}_{-3.39}^{+3.26}$	$-{4.89}_{-3.54}^{+3.43}$	$-{3.43}_{-4.41}^{+2.50}$	$-{2.26}_{-3.88}^{+1.71}$
$\mathrm{log}\,{{\rm{P}}}_{0}$	0.0	$-{0.16}_{-0.80}^{+1.32}$	$-{0.19}_{-0.71}^{+1.03}$	${0.03}_{-0.74}^{+1.16}$	${0.45}_{-0.85}^{+1.01}$
${R}_{{\rm{p}}}$	1.0	${1.13}_{-0.50}^{+1.60}$	${1.13}_{-0.48}^{+1.27}$	${1.02}_{-0.38}^{+1.10}$	${0.80}_{-0.19}^{+0.81}$
$\mathrm{log}\,g$	0.99	${1.42}_{-0.82}^{+0.42}$	${1.45}_{-0.75}^{+0.41}$	${1.41}_{-0.83}^{+0.43}$	${1.26}_{-0.83}^{+0.52}$
$\mathrm{log}\,{A}_{s}$	−1.3	$-{0.89}_{-0.74}^{+0.63}$	$-{0.84}_{-0.70}^{+0.56}$	$-{0.95}_{-0.68}^{+0.64}$	$-{0.76}_{-0.79}^{+0.53}$
$\mathrm{log}\,{p}_{t}$	−0.22	$-{1.26}_{-0.52}^{+0.98}$	$-{1.24}_{-0.55}^{+0.83}$	$-{1.23}_{-0.57}^{+0.84}$	$-{0.84}_{-0.73}^{+0.89}$
$\mathrm{log}\,{dp}$	−1.0	$-{1.75}_{-0.87}^{+1.16}$	$-{1.70}_{-0.87}^{+1.12}$	$-{1.46}_{-1.06}^{+1.14}$	$-{1.49}_{-1.03}^{+1.49}$
$\mathrm{log}\,\tau$	1.0	${0.03}_{-1.39}^{+1.33}$	${0.05}_{-1.40}^{+1.42}$	${0.61}_{-1.55}^{+0.94}$	${0.99}_{-1.44}^{+0.73}$
$\mathrm{log}\,{f}_{c}$	−0.3	$-{1.41}_{-1.07}^{+0.96}$	$-{1.42}_{-1.07}^{+1.00}$	$-{0.82}_{-1.28}^{+0.60}$	$-{0.58}_{-0.97}^{+0.41}$

Download table as:  ASCIITypeset image

For R = 70 at S/N = 5, Figure 11 shows that there is only a weak detection of P₀ and a detection of R_p, which merely suggests that the planet has an atmosphere and is not a giant planet. As S/N increases to 10, the O₃ posterior distribution has a weak peak near the fiducial value, and the gas is only weakly detected. Once the S/N is equal to 15, we weakly detect ${{\rm{H}}}_{2}{\rm{O}}$, O₃, and O₂. At an S/N of 20, it is possible to detect each of ${{\rm{H}}}_{2}{\rm{O}}$, O₃, and O₂. At this S/N, the oxygen mixing ratio is estimated to be above roughly 10⁻³, indicating that we are unable to determine if O₂ is a major atmospheric constituent (i.e., present at the 1% level or more). Gravity (and, thus, planetary mass) remains undetected at all S/Ns, similar to the findings of Lupu et al. (2016). The surface albedo is unconstrained (or worse) at all S/Ns but shows a weak bias toward a higher value of A_s ≈ 0.3 ($\mathrm{log}{A}_{{\rm{s}}}\approx -0.5$) at the highest S/Ns, which is likely due to the relatively large error bars at red wavelengths (driven primarily by low detector quantum efficiency) where we have the most sensitivity to the surface. We are able to get weak detections of τ and f_c, which are shown in Figure 10 to be correlated. Yet with these posteriors, we cannot rule out scenarios without cloud cover. We note that the drop-off in the posteriors of p_t and dp at higher pressure values likely result from the limiting conditions that the cloud base cannot extend below the surface pressure and the upper limit of the P₀ prior. The improved S/N leads to a posterior more concentrated around the true value for p_t, dp, and R_p. For improved constraints on cloud properties, it may be beneficial to observe time variability with photometry (e.g., Ford et al. 2001) or use polarimetry (e.g., Rossi & Stam 2017).

At a higher spectral resolution (R = 140), the improvement in detections and constraints begins at a lower S/N, as illustrated by Figure 12. Gravity remains undetected for all S/Ns. At an S/N equal to 5, P₀ and R_p have a weak detection and a detection, respectively. At S/N = 10, it is possible to detect ${{\rm{H}}}_{2}{\rm{O}}$, O₃, and O₂. As with the R = 70 case, the surface albedo is unconstrained (or worse) at all S/Ns, and, at the highest S/Ns, the model is biased toward A_s ≈ 0.3 (as with R = 70). Moving to S/N = 15 adds a constraint to R_p, P₀, and O₃, as well as weak detections of cloud parameters. Increasing the S/N to 20 does not dramatically change the posterior distributions, although the posteriors for ${{\rm{H}}}_{2}{\rm{O}}$ and O₂ become narrow enough to offer constraints. Here the constraint on O₂ suggests it is a major constituent in the atmosphere. In spite of the generous S/N, though, the 1σ uncertainties on the gas mixing ratios are not more precise than roughly an order of magnitude (see Table 5).

Considering both R = 140 and R = 70, we see that S/N = 5 data offer very little information about the planetary atmosphere. In the case of R = 140, S/N = 10 data offer detections but no constraints, and S/N = 20 data are required to constrain all included gas species. In other words, the conclusions we would draw about the planet (e.g., the amount of gases, the bulk and cloud properties) improve significantly between S/N = 10 and 20. With R = 70, the boost from S/N = 10 to 15 provides weak detections of key atmospheric and surface parameters, and S/N = 20 data offer detections but few constraints (i.e., except on planetary radius).

For the WFIRST rendezvous data sets, we are able to infer very little information at an S/N of 5 or 10 except for weak detection of surface pressure and a detection of the planetary radius. All gases remain undetected at these S/Ns. The posterior distributions for most parameters do not vary much as S/N improves, although there are weak detections of cloud optical depth and fractional coverage at the highest S/Ns. Like all previous cases, we do not detect the surface gravity. At S/N = 15 and 20, the detection of f_c is unable to rule out scenarios with little cloud cover. To obtain weak detections of the atmospheric gases, we require an S/N of 20, but even here, the posteriors have tails that extend to near-zero mixing ratios.

To compare the performance of a WFIRST rendezvous scenario against HabEx/LUVOIR scenarios at R = 70 and 140, we plot the posterior distributions of the parameters for the S/N = 10 results from the WFIRST rendezvous and R = 70 and 140 in Figure 15. While this comparison sheds light on the corresponding trade-off in terms of parameter estimation for the same S/N, these cases do not represent equal integration times, which scale with resolution and S/N. If the dominant noise source does not depend on resolution (e.g., detector noise), the cases of R = 140 at S/N = 10, R = 70 at S/N = 20, and R = 50 at S/N = 28 would be roughly equal in integration time. However, if the dominant noise source does depend on resolution (e.g., exozodiacal dust), the cases of R = 140 at S/N = 10, R = 70 at S/N = 14, and R = 50 at S/N = 17 would have roughly equivalent integration times. Tables 7–9 allow approximate comparisons of these different scenarios, excluding a WFIRST rendezvous scenario at high S/N = 28 that we have not considered.

Zoom In Zoom Out Reset image size

From Figure 15, we see that the performance of the WFIRST rendezvous retrieval is similar to that of R = 70 at S/N = 10. The noticeable difference is a weak detection of O₃ with R = 70. Because we adopt the photometric setup from WFIRST Design Cycle 7 through the shorter wavelengths, the data do not provide complete spectroscopic coverage across the significant O₃ feature from 0.5 to 0.7 μm, as in the case of HabEx/LUVOIR simulated data. Figure 5 shows the sampling of the forward-model spectrum for the three types of data sets we considered. We compare the spectral fits in Figure 14 and note the much wider spread in the possible fits for wavelengths shorter than 0.6 μm for the WFIRST rendezvous versus R = 70 or 140, which have continuous coverage in the full range. The R = 140, S/N = 10 data set was able to offer detections of all atmospheric gases, setting it apart from the other two. We stress, however, that constraints were only found at S/N = 20 and R = 140.

Our parameter estimations are likely to be optimistic as a consequence of our adoption of nonrandomized spectral data points in our faux observations. Thus, the requisite S/Ns for detection detailed above should be seen as lower limits. Ultimately, our decision to use nonrandomized data points stemmed from computational limitations (preventing us from running large numbers of randomized faux observations for each of our R–S/N pairs) and a desire to avoid the biases that can occur from attempting to make inferences from retrievals performed on a single randomized faux observation (Lupu et al. 2016).

However, we deemed it necessary to investigate the consistency of our findings with respect to different noise instances. To work within our computational restrictions, we realized that cases such as R = 70 with S/N = 5 yielded little detection information for any parameter, even in the ideal scenario of nonrandomized data. We then decided to focus on two "threshold" cases based on the results from the nonrandomized data: R = 140 with S/N = 10 and R = 70 with S/N = 15. We ran 10 noise instances of these two cases where it is likely the optimistic nonrandomized data make the difference between detection and constraint for several parameters (see Tables 7 and 8).

Each noise instance is run for at least 10,000 steps in emcee. Figure 16 shows all the individual posteriors for the gas mixing ratios from each noise instance for R = 70, S/N = 15. We highlight the posteriors from one "outlier" case where there is no oxygen detection. The corresponding set of data points is shown as well. This highlights the fact that single noise instances can mislead our interpretation and the benefit of having many noise instances run to obtain a more comprehensive understanding of the state of an atmosphere.

Zoom In Zoom Out Reset image size

To summarize the noise instance results, we concatenate samples from the last 1000 steps in each noise instance and construct an averaged set of posteriors. We are able to do this because the noise instances are equally likely, having been drawn in the same manner from a Gaussian with set parameters (i.e., the same S/N as the standard deviation). In Figure 17, we plot the combined posteriors of the 10 noise instances of R = 70, S/N = 15 and compare them to the posterior from the last 5000 steps of the nonrandomized data case. We illustrate the same comparison for R = 140, S/N = 10 in Figure 18. We overplot the truth values, as well as the 68% confidence interval and median value, for each parameter from the combined noise-instances posterior and the nonrandomized data posterior.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For all parameters in both the R = 70 and R = 140 cases, we find that the average posterior from the 10 noise instances agrees with the posterior from the nonrandomized data set qualitatively. Their medians and 68% confidence interval ranges are also similar with significant overlap. The overall conclusions we can draw from the average posteriors do not appear to differ much from those using the nonrandomized data set posteriors.

Future space-based direct-imaging missions will have a diversity of goals for exoplanet studies and will likely emphasize the detection and characterization of Earth-like exoplanets. For the detection of oxygen and ozone—which are key biosignature gases—in the atmospheres of Earth twins, our results indicate that spectra at a minimum characteristic S/N of 10 will suffice if at R = 140, while data at an S/N of at least 15–20 would be needed at R = 70. For a WFIRST rendezvous-like observing setup, these gases would only be weakly detected, even at an S/N of 20. Methane, which is another important biosignature gas, has no strong signatures in the visible wavelength range for the modern Earth, so we did not consider detection of this gas. Thus, we could not use our simulated data and retrievals to argue for detections of atmospheric chemical disequilibrium (Sagan et al. 1993; Krissansen-Totton et al. 2016).

Key habitability indicators include atmospheric water vapor and surface pressure. Detecting the former requires an S/N of 15–20 at R = 70 but only an S/N of 10 at R = 140. Surface pressure can be constrained to within an order of magnitude for S/N ≳ 15 at R = 140, although the overall lack of temperature information in these reflected-light spectra would make it impossible to use pressure/temperature data to argue for habitability (Robinson 2017). Surface temperature information may then need to come from climate-modeling investigations that are constrained by retrieved gas mixing ratios.

For all of our observing setups, the data yield detections of, and in some cases constraints on, the planetary radius. Except at an S/N of 20 for R = 70 or S/N > 15 for R = 140, the posterior distributions are not well enough constrained to distinguish an Earth/super-Earth (R_p < 1.5 R_⊕) from a mini-Neptune based on size alone, although the data do rule out planetary sizes larger than Neptune. Additional atmospheric information (e.g., composition) could potentially be used to help distinguish between terrestrial planets and mini-Neptunes. These findings are consistent with the gas giant–focused work of Nayak et al. (2017), who noted that observations at multiple phase angles can also help to better constrain planetary size. Our overall lack of surface gravity constraints, paired with the weak constraints on planet size, implies that we do not have a constraint on the planetary mass. Follow-up (or precursor) radial velocity observations (or, potentially, astrometric observations) could offer additional constraints on planet mass.

We can make rough comparisons of our R–S/N results to those of Brandt & Spiegel (2014), who used minimally parametric models to investigate detections of O₂ and ${{\rm{H}}}_{2}{\rm{O}}$ for Earth twins. These comparisons are not direct, however, as Brandt & Spiegel (2014) were fitting for fewer parameters (8 versus our 11) and also only assumed that the S/N was proportional to planetary reflectance (versus our more complicated scaling, as shown in Figure 6). For O₂, Brandt & Spiegel (2014) found R = 150 and S/N = 6 for a 90% detection probability, which is consistent with our R = 140 posteriors moving from a nondetection at S/N = 5 to a detection at S/N = 10. When investigating ${{\rm{H}}}_{2}{\rm{O}}$, Brandt & Spiegel (2014) found R = 40 and S/N = 7.5 or R = 150 and S/N = 3.3 for a 90% detection probability. Using Figure 6 to scale our S/Ns to 890 nm (i.e., to the continuum just shortward of the 950 nm water vapor band), at R = 50, we only find a weak detection of ${{\rm{H}}}_{2}{\rm{O}}$ for S/N_890nm  = 10, and at R = 140, we transition from a water vapor nondetection to a detection between an S/N_890nm of 2.5–5. Taken all together, these comparisons indicate that we agree with Brandt & Spiegel (2014) at higher spectral resolution (R = 140–150) but that detection of ${{\rm{H}}}_{2}{\rm{O}}$ at lower spectral resolution (R = 50) will likely require higher S/Ns than originally indicated.

The discussion above emphasizes mere detections, not constraints (which, again, we define as having peaked posterior distributions with 1σ widths less than an order of magnitude). While uncertain, we anticipate that characterization of climate, habitability, and life likely require constraints, not simple detections. Here, as is shown in Table 8, only R = 140 and S/N = 20 observations offer the appropriate constraints. Thus, future space-based high-contrast imaging missions with goals of characterizing Earth-like planetary environments are likely to need to achieve R = 140 and S/N = 20 observations (or better). Of course, combining near-infrared capabilities, which would provide access to additional gas absorption bands, may help loosen these requirements.

Several key assumptions adopted in this study warrant further comment. First, as noted earlier, we do not retrieve on planetary phase angle and planet–star distance, both of which influence the planet-to-star flux ratio. Thus, in effect, we are assuming that the planetary system has been revisited multiple times for photometric and astrometric measurements, such that the planetary orbit is reasonably well-constrained (i.e., that the orbital distance and phase angle are not the dominant sources of uncertainty when interpreting the observed planet-to-star flux ratio spectrum). If the orbit is not well-constrained, Nayak et al. (2017) showed that strong correlations can exist between the retrieved phase angle and the planet radius.

Second, we have assumed detector-dominated noise and a quantum efficiency appropriate for the WFIRST/CGI for all of our observational setups. While this is likely a fair assumption for our WFIRST rendezvous studies, it is likely that detector development will lead to major improvements in instrumentation for a HabEx/LUVOIR-like mission. Here the rapid decrease into the red due to detector quantum efficiency may not be as dramatic, implying that spectra would have relatively more information content at red wavelengths as compared to the present study. Furthermore, a HabEx/LUVOIR-like mission may no longer be in the detector-dominated noise regime. In the limit of noise dominated by astrophysical sources (e.g., exozodiacal light or stellar leakage), the S/N only varies as $\sqrt{q{ \mathcal T }{B}_{\lambda }}$.

Finally, we adopt a relatively simple parameterization of cloud 3D structure. Specifically, we allow for only a single cloud deck in the atmosphere, and we then permit these clouds to have some fractional coverage over the entire planet. This parameterization of fractional cloudiness implies uniform latitudinal and longitudinal distribution of patchy clouds. In reality, clouds on Earth have a complex distribution in altitude, latitude, and longitude (Stubenrauch et al. 2013), and variations in time also have an observational impact (Cowan et al. 2009; Cowan & Fujii 2017). However, given the overall inability of our retrievals to constrain cloud parameters (at least at the S/Ns investigated here; see also Lupu et al. 2016; Nayak et al. 2017), it seems challenging for future space-based exoplanet characterization missions to detect (or constrain) more complex cloud distributions with the types of observations studied here and data of similar quality.

Our current forward model is able to include both CO₂ and CH₄, although we did not retrieve on these gases in the current study due to their overall lack of strong features in the visible wavelength range for modern Earth. However, these species do have stronger features in the near-infrared wavelength range. As both of the HabEx and LUVOIR concepts are considering near-infrared capabilities, it will be essential to extend our current studies to longer wavelengths and investigate whether or not constraints on additional gases (i.e., beyond water, oxygen, and ozone) can be achieved at these wavelengths.

Additionally, given the likely huge diversity of exoplanets that will be discovered by future missions (and that have already been identified and studied by Kepler, Hubble, and Spitzer), it will be necessary to extend our parameter estimation studies to include a wider range of worlds. Both super-Earths and mini-Neptunes are more favorable targets for a WFIRST rendezvous mission and may also be easier targets for HabEx/LUVOIR-like missions. Our forward model is already capable of simulating these types of worlds, and we anticipate emphasizing a variety of exoplanet types in future studies. Such future studies may also include retrievals on planetary phase angle, which would be relevant to observing scenarios where the planetary orbit is poorly constrained.