Exoplanet Classification and Yield Estimates for Direct Imaging Missions
Ravi kumar Kopparapu Eric Hebrard Rus Belikov Natalie M. Batalha Gijs D. Mulders Chris Stark Dillon Teal Shawn Domagal-Goldman Avi Mandell




1NASA Goddard Space Flight Center, 8800 Greenbelt Road, Mail Stop 699.0, Building 34, Greenbelt, MD 20771
2Department of Astronomy, University of Maryland College Park, College Park, MD
3NASA Astrobiology Institute's Virtual Planetary Laboratory, P.O. Box 351580, Seattle, WA 98195, USA
4Blue Marble Space Institute of Science, 1001 4th Ave, Suite 3201, Seattle, Washington 98154, USA
5School of Physics and Astronomy, University of Exeter, EX4 4QL, Exeter, UK
6NASA Ames Research Center
7Lunar and Planetary Laboratory, University of Arizona, Tucson, AZ 85721, USA
8Earths in Other Solar Systems Team, NASA Nexus for Exoplanet System Science, USA
9Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218 





Future NASA concept missions that are currently under study, like Habitable Exoplanet Imaging Mission (HabEx) & 
Large Ultra-Violet Optical Infra Red (LUVOIR) Surveyor, would discover a large diversity of exoplanets.
We propose here a classification scheme that distinguishes exoplanets into different categories based on their size and incident 
stellar flux, for the purpose of providing the expected number of exoplanets observed (yield) with direct imaging missions. 
The 
boundaries of this classification can be computed using the known chemical behavior of gases and condensates at different pressures
and temperatures in a planetary atmosphere. In this study, we initially focus on condensation curves for sphalerite ZnS, 2o, 
2 and 4. The order in which these species condense in a planetary atmosphere define the boundaries between different
classes of planets.  Broadly, the planets are divided into rocky (0.5 - 1.0R_⊕), super-Earths 
(1.0 - 1.75R_⊕), sub-Neptunes (1.75 - 3.5R_⊕), sub-Jovians (3.5 - 6.0R_⊕) and 
Jovians (6 - 14.3R_⊕) based on their planet sizes, and 'hot', 'warm' and 'cold' based on the incident stellar flux. We
then calculate planet occurrence rates within these boundaries for different kinds of exoplanets, η_planet, 
using the community co-ordinated results of NASA's Exoplanet 
Program Analysis Group's Science Analysis Group-13 (SAG-13). These occurrence rate estimates are in turn used to estimate the 
expected exoplanet yields for direct imaging missions of different telescope diameter. 



Exoplanet Classification and Yield Estimates for Direct Imaging Missions
    Ravi kumar Kopparapu1,2,3,4, Eric Hébrard1,5, Rus Belikov6,  Natalie M. Batalha6, Gijs D. Mulders7,8, Chris Stark9, Dillon Teal1, Shawn Domagal-Goldman1,3, Avi Mandell1
    December 30, 2023
============================================================================================================================================================================






§ INTRODUCTION

The discoveries of exoplanets over the last two decades has revealed planetary bodies of various sizes and masses around other stars <cit.>. More specifically, the location of these exoplanets around their host star has a significant influence not only on the prospects of their detectability, but also on the atmospheric chemical composition and the capability of characterizing these atmospheres. Furthermore, the relatively large sample of nearly 3400 confirmed planets and another ∼ 4700 planet candidates to-date[http://exoplanetarchive.ipac.caltech.edu/] has enabled us to calculate exoplanet occurrence rates in our galaxy <cit.>. These initial estimates are dominated by close-in planets due to the sensitivity of the detection techniques and search pipelines. Nevertheless, these studies made a crucial and a significant leap in understanding planet diversity, and paved a way for comparative planetology of exoplanets. 

Several of the above mentioned studies have also focused on obtaining an estimate of the fraction of stars that have at least one terrestrial mass/size planet in the habitable zone (HZ), or η_⊕. 
Estimates of η_⊕ for Sun-like stars have been calculated by the data collected from the Kepler mission. Earlier estimates ranged from 0.02 <cit.> to 0.22 <cit.> for GK dwarfs, but more recent analyses <cit.> imply that systematic errors dominate.  
For M-dwarfs, η_⊕ is estimated to be ∼ 20% 
on an average <cit.>. Apart from the general curiosity of finding how common are Earth-like planets in our galaxy, the focus on η_⊕  has a more practical application: It can be used in the design of direct imaging missions, like the concept studies under consideration HabEx[http://www.jpl.nasa.gov/habex/], the Habitable Exoplanet Explorer; and LUVOIR[https://asd.gsfc.nasa.gov/luvoir/], the Large UV-Optical-InfraRed surveyor
 with the goal of detecting biosignatures, and also in calculating `exo-Earth candidate yield', the number of potentially habitable extrasolar planets (exoEarth candidates) that can be detected and spectroscopically characterized <cit.>. 


Crucial to these estimates are the location of the main-sequence HZ, which have been studied by both 1-D and 3-D climate models <cit.>. With some exceptions for certain types of planets <cit.>, there has not been an overarching way to classify planets beyond the HZ. The lack of a systematic way to classify exoplanets in general, combined with the allure of planets within it, has led to direct imaging mission yield analyses that focus on HZ planets to the exclusion of everything else <cit.>. While some mission studies have attempted to classify the non-HZ planets into hot, warm and cold planets that a mission would discover[https://exoplanets.nasa.gov/exep/studies/probe-scale-stdt/], the boundaries for the classification are arbitrarily fixed without giving consideration to chemical behavior of gases and condensates in a planetary atmosphere.  Classifying different size planets  based on the  transition/condensation of different species 
<cit.>
provides a physical motivation in estimating exoplanet mission yields, separate from exo-Earth candidate yields. 

In the search for exo-Earth candidates, we will undoubtedly detect a multitude of brighter planets. According to <cit.>, for an 8m size telescope, the number of exo-Earth candidates detected is ∼ 20 (see Fig. 4 in Stark et al. 2014), although this is strongly dependent on the value of η_⊕. 
At the same time, the number of stars observed to detect these exo-Earth candidates is ∼ 500. If we assume that, on an average, every star has a planet of some 
size <cit.>, then there are  ∼ 500 exoplanets of all sizes that can be observed. Not considering the ∼ 20 exo-Earth candidates, the bulk of the exo-planets will fall into `non-Earth' classification, without any distinguishing features between them. This provides an additional motivation to devise a scheme based on planetary size and corresponding atmospheric characteristics of exoplanets.   

Some work has been done at the theoretical level 
to derive the radiative response of an irradiated atmosphere
<cit.>, but these analytical tools have not been either used to derive 
any general boundaries nor tied to planet occurence estimates,
 and were not designed with that intent in mind. This highlights the need for a theory-based system to classify planets beyond the habitable zone for the purposes of understanding the diversity of worlds future missions could explore. Such a system should be based on the properties we can measure today, primarily size and orbital information, and the boundaries in the scheme should divide planets with major differences in the properties that would be observable with current and future missions.


Fortunately, much of the theory needed for such a classification scheme already exists. There has been significant progress in understanding how the size of a planet is a major control on composition, and therefore on future observables <cit.>. This includes work on the relationship between size, density, and bulk composition <cit.>.  


While the exact values for the boundaries of mass-radius change based on the specific analysis or theoretical technique, 
there is growing evidence of structure in the occurrence rate distribution that suggests compositional aggregates:
 1) small, rocky worlds whose bulk composition and behavior is dominated by Fe, Mg and Si species; 2) planets with
 Fe/Mg/Si cores but significant gas envelopes consisting of H/He, 4, 3 ices; and 3) gas giants whose bulk composition and behavior is dictated almost exclusively by its volatiles.

Similarly, there have been significant advances in our understanding of how the orbital separation of non-HZ planets could affect the chemical composition of the atmosphere <cit.>. A constant theme across these studies is the influence of clouds. As a planet moves further from its host star, its atmosphere will cool and lead to condensation of progressively less volatile chemicals in the atmosphere. This condensation would create a cold trap and an associated cloud deck. The result of this is a significant change in spectral properties: as the condensing species would be trapped at or below the cloud deck, the cloud deck itself would absorb and scatter light, causing preferential sampling of the layers at or above the cloud deck. Multiple "onion-like" cloud decks can form as sequentially less volatile species condense at higher altitudes for planets with greater star-planet separation distances and correspondingly lower levels of incoming stellar flux. This process can be observed in detail in the gas giants of our own Solar System <cit.>, and has been postulated to be a driver for the atmospheric structure and observable properties of exoplanet atmospheres <cit.>. We also note that the habitable zone itself has been defined in a manner consistent with this, as the instellations (stellar incident flux) at which liquid water clouds form but carbon dioxide clouds do not <cit.>.

This prior work enables the construction of an overarching scheme for identifying classes of planets. This scheme could apply to all worlds, regardless of whether they are rocky or gaseous (or something in between). And it would be both based on the current observable properties, and prediction of major transitions in future observables. In short, this represents a comprehensive means of predicting the diversity yields of future planet characterization missions. Below we discuss in more detail how we simulate the processes underlying this classification scheme. These simulations define the boundaries between different planet classes, for which we calculate occurrence rates based on prior exoplanet detection missions. The occurrence rates allow us to simulate exoplanet yields - not just for HZ planets but for a diversity of different kinds of worlds. Finally, we close with a discussion of the caveats of this approach, and the implications of this scheme for future missions.

 

§ A NEW CLASSIFICATION SCHEME

        


A planet size - and the relationship between its size and mass - appears to be primarily driven by volatile inventory. For example, the atmospheric composition of larger planets is predominantly  H_2/He, while smaller planets can have a mixture of 4, 2, 2o and 3. High-temperature atmospheres, such as hot Jupiters, should have their chemistry - and therefore their spectral features - determined primarily by equilibrium chemistry.  Low-temperature atmospheres will have chemistry dictated by photochemistry, but this will be secondary to determining what species are condensing in their atmospheres. The exception to this - which we will discuss later - is for photochemical aerosols, which could have a major impact in the same manner that clouds do.


The chemical behavior of gases and condensates in a planetary atmosphere can be determined as a function of pressure, temperature, and metallicity. Using results from Lodders and Fegley (2002) and Visscher et al. (2006) adapted with solar abundances taken from Lodders (2010), we have computed the condensation curves for sphalerite ZnS, 2o, 2 and 4 as a function of pressure and temperature for systems with a solar metallicity. 
Pressure-temperature profiles of planetary atmospheres are tightly related to the incoming stellar flux. We define the boundaries between our different selected planetary cases as the stellar fluxes for which these four species condense out.
For instance, ZnS clouds have been considered as possible condensates in hot exoplanet atmospheres <cit.>, so the location (or the stellar flux) at which ZnS clouds form in a planetary atmosphere denotes the first boundary for our planet classification. Moving further away from the star, at relatively lower stellar fluxes, 2o starts condensing in the atmosphere. This, then, becomes the next boundary for grouping different planets. Between these two boundaries is where one would expect to find ZnS mineral clouds and 2o in a gaseous state. Continuing to lower fluxes, 2 and 4 condensates bracket final boundaries. The results are independent of any particular model atmospheres, and in principle, any pressure-temperature profile may be superimposed on these condensation curves to find the equilibrium composition along the profile. In particular, the intersection of a particular pressure-temperature profile with one of these condensation curves indicates the pressure and temperature at which the respective species condenses out in the planetary atmosphere considered. We have investigated the different incident stellar fluxes, or instellations, for which condensates could form within the regions of an exoplanet system that will be probed by future direct imaging missions. Specifically, we simulate the star-planet separations for which ZnS, 2o, 2 and 4 would condense out in planetary atmospheres. Other metallic clouds can condense at distances closer to the star than the ZnS condensation line, and other volatiles (e.g., NH_3) can condense at orbiting distances beyond the 4 condensation line, but such worlds are likely undetectable by future direct imaging missions so are not simulated here.

 We have considered six different planetary size boundaries: 0.5 R_⊕, 1.0 R_⊕, 1.75 R_⊕, 3.5 R_⊕, 6.0 R_⊕, and 14.3 R_⊕ in our grid. 
These boundaries represent, respectively: the radius (0.5 R_⊕) at which planets in the habitable zone appear to not have a sufficient gravity well to retain atmospheres <cit.>;
the “super-Earths” (1 - 1.75 R_⊕) and “sub-Neptunes” (1.75 - 3.5 R_⊕), as defined by
<cit.> (see section 4.4) based on the observed gap in the radius distribution of small planets 
with orbital periods shorter than 100 days;
the assumed upper limit on Neptune-size planets (6R_⊕) based on the small peak in the radius distribution
from <cit.>;


and the radius past which planets transition to brown dwarf stars <cit.>. We have computed the 
corresponding pressure-temperature atmospheric profiles using (1) the non-grey analytical model of 
Parmentier and Guillot (2014) with the coefficients from Parmentier et al. (2015) and the Rosseland opacity functional 
form of Valencia et al. (2013), and (2) the grey analytical model of Robinson and Catling (2012, 2014), both modified 
to take the planetary size and instellation as unique input parameters. We have used the 
Robinson and Catling (2012, 2014) model for planets with radius smaller than 3.5 R_⊕ at low instellations, 
and the Parmentier and Guillot (2014) model for planets with a radius of 14.3 R_⊕ at high instellations. We 
have assumed an internal temperature T_int=0 and mass-radius relations taken from Weiss and Marcy (2014) for 
planets with radius smaller than 3.5 R_⊕. We have assumed an internal temperature T_int=100 and the 
density calculated from Mass-radius relation from <cit.> for planets with radius equal to 
14.3 R_⊕. Assuming a planetary mass of 0.414M_J (131M_⊕), the density is ∼ 0.25
g.cm^-3.  We have 
considered L = 0.95 L_⊙ and M = 0.965M_⊙ for the parent star[The values of the stellar
luminosity and mass are obtained as follows: We downloaded the confirmed and candidate Kepler catalog from NEXSCI,
 found the median values of stellar luminosity and mass for each data set, and then took the average value of 
luminosity and mass from these median values.
Median Luminosity for candidate planet list: 1.057;
Median Stellar mass for candidate planet list: 0.97;
Median Luminosity for confirmed planet list: 0.86;
Median Stellar mass for confirmed planet list: 0.96;
Average luminosity of confirmed & candidate: (1.057 + 0.86)/2 = 0.95;
Average stellar mass of confirmed & candidate: (0.97+0.96)/2 = 0.965].

As an illustrative example, we show in Fig. <ref> the condensation curves for the four different species we 
considered (solid black), along with temperature profiles for two different size planets: 0.5R_⊕ and 
14.3R_⊕, at different incident stellar fluxes (`instellation'). The figure shows that ZnS would condense 
out at ∼10mb in the atmosphere of a highly irradiated (220 times Earth flux) 14.3 Earth radius planet, a typical 
hot-Jupiter, whereas 4 would condense out in the atmosphere of a 0.5  Earth radius planet receiving only 
∼ 1/280^th the flux Earth is receiving. Following this procedure, we have derived the radius and stellar flux 
where other gaseous species condense in the atmosphere. Table 1 provides the corresponding data of the planetary 
radius and stellar flux boundaries that can be used for classifying planets into different regimes.  These
boundaries are parameterized in Eq.(<ref>) and are also available as an online calculator at: 
<http://www3.geosc.psu.edu/ ruk15/extrasolar/>.

 [h!]

 Planetary radii and stellar flux values at which the given species condense in a planetary atmosphere, 
for a star with L = 0.95 L_⊙ and M = 0.965M_⊙. These 
limits form the boundaries for classifying planets into different categories to calculate exoplanet yield estimates. See Fig. 2 and Section 2.2.



        Stellar Flux (Earth flux)         

  Radius (R_⊕)     ZnS      2o     2     4
 
 0.5   182   1.0   0.28   0.0035

 1.0   187   1.12   0.30   0.0030

 1.75   188   1.15   0.32   0.0030

 3.5   220   1.65   0.45   0.0030

 6.0   220   1.65   0.40   0.0025

 14.3   220   1.7   0.45   0.0025

 





    F(R_p)_i   =    a_i x^5 + b_i x^4 + c_i x^3 + d_i x^2 + e_i x + f_i


where F(R_p)_i is the stellar flux, normalized to the current Earth flux (1360 Wm^-2), at which species
i = (ZnS, 2o, 2, 4) condenses on a planet with radius R_p, and x = R_p/R_⊕ given in
Table <ref>. The 
coefficients for the four condensing species are given in Table <ref>.

[h!]

 Coefficients to be used in Eq.(<ref>).



    a   b   c   d   e   f  

  ZnS     0.0010338041    - 0.0255230451     0.1858822989     - 0.4990171468     0.5690844110     1.6385396777 

 2o     0.0017802416     - 0.0443704003     0.3302966853     - 0.9299682851     1.1366785108     0.6255832476 

 2     0.0002947546     - 0.0070583509     0.0483928147     - 0.1198359100     0.1477297602     0.2304769313 

 4     -0.0000033096     0.0000889715     -0.0007644190     0.0026719183     -0.0038305535     0.0048373922 

 









 §.§ Application of the Classification Scheme to Obtain Planet Occurrence Rates



Extending the insights obtained from Fig. <ref> and Table 1, it is then possible to define various ‘zones’ as a function of stellar flux and planetary radius. This should qualitatively apply across all stellar types and the entire field of planets, even if the quantitative positions of the boundaries change due to, for example, the age of the system and the amount of internal heat released from planets, or the stellar energy distribution (SED) and its relationship to planetary albedo <cit.>. For a full consideration of such caveats, see section 3.1. This particular framing of the parameter space will allow us to calculate the occurrence of different kinds of planets within each zone based on their condensation conditions, as both radius and flux are measurable quantities.

We have performed such a calculation of η_planet, the fraction of stars that have a planet within one of the zones defined by a condensing species. 
As an illustrative example of how these boundaries can be used to calculate the occurrence of planets, we apply these criteria to the preliminary parametric model introduced by one of NASA's Exoplanet Program Analysis Group (ExoPAG) science analysis studies (SAG13). A detailed discussion of the SAG13 model is outside the scope of this paper, but we will summarize the most critical points.

The SAG13 model is based on a simple meta-analysis of planet occurrence rates from many different individual publications and groups. Specifically, the SAG13 group collected tables of occurrence rates calculated over a standard grid of 
planet radius, period, and stellar type. A full description of the grid is as follows.
 The i^th bin in the planet radius is defined as the interval:


    R_i   =    [1.5^i-2, 1.5^i-1) R_⊕

This implies the following bin edges: [0.67, 1.0, 1.5, 2.3, 3.4, 5.1, 7.6, 11, 17,...] R_⊕.

The j^th bin in the planet period is defined as:

    P_j   =    10 . [2^j-1, 2^j) days

This implies the following bin edges:
[10, 20, 40, 80, 160, 320, 640,...] days


Data and models from peer-reviewed publications (Petigura et al. 2013, ForemanMackey et al. 2014, Burke et al. 2015, 
Traub 2015, Dressing & Charbonneau 2013) were integrated over the standard grid, and supplemented by several 
unpublished tables from the 2015 Kepler “hack week” which were based on Q1-Q17 DR24 catalog, Kepler completeness 
curves, and data products at the time.

 However, for our current work, we did not use the SAG13 standard grid mentioned above because the SAG13 grid does not 
represent the condensation sequences described in the previous section. Instead, we took the stellar fluxes from Table 1 
where species condensation happens, and calculated the corresponding orbital periods based on the stellar mass 
(0.965) M_⊙ and luminosity (0.95) L_⊙ described in footnote 5. 
It should be noted that the SAG-13 grids are available for different spectral types.  Herein, our work focuses on the G dwarf population and employs the corresponding grids.


 The SAG13 submissions were then processed as follows. First, within each spectral type, the sample geometric mean 
(μ_i,j) and variance (σ_i,j^2) was computed in each (i,j)-th bin of the period-radius grid, across the 
different submissions. The mean values μ_i,j formed a “baseline” table of occurrence rates. “Optimistic” and 
“pessimistic” tables were also defined by using the μ_i,j± 1 σ_i,j values for each (i,j)-th bin.

SAG13 then fit a piecewise power law to the “pessimistic”, “baseline”, and “optimistic” combined tables. 
The power law had the following form:



    ∂^2N(R,P)/∂ lnR ∂ lnP   =   Γ_iR^α_iP^β_i


  For optimistic case, Γ_i = [1.06, 0.78], α_i = [-0.68, -0.82], β_i = [0.32, 0.67]; for pessimistic case
Γ_i = [0.138, 0.72], α_i = [0.277, -1.56], β_i = [0.204, 0.51]; for baseline, Γ_i = [0.38, 0.73], α_i = [-0.19, -1.18], β_i = [0.26, 0.59]

 The break between two pieces of the power law was set at 3.4 R_⊕ (following Burke et al. 2015), hence the two values
for the coefficients, and a least squares fit 
was performed separately to each of the pieces. Similarly to the mean and variance above, logarithms of occurrence rates were used 
when performing the least squares of log occurrence rates, rather than actual occurrence rates, in order to properly balance the 
effects of small and large occurrence rates. This resulted in “pessimistic”, “baseline”, and “optimistic” parametric models. 
These models were then integrated across the planet parameter boundaries described in this paper.

It should be stressed that community-sourced data do not represent independent measurements or estimates of scientific quantities, 
so that the SAG13 sample mean and variances should not be interpreted as a formal mean and uncertainty of exoplanet occurrence rates. 
Rather, they simply represent one possible way to measure the state of knowledge as well as the disagreement on the rates within the 
occurrence rate community. In other words, the SAG13 “pessimistic”, “baseline”, and “optimistic” cases refer to the typical 
pessimistic, average, and optimistic submissions within the SAG13 community survey, rather than formal scientific results. 

Alternative ways of combining SAG13 results are also possible, such as: including only peer-reviewed submissions, including 
submissions based only on different catalogs, removing outliers, etc. A detailed analysis of this is beyond the scope of this paper, 
but as a general rule, combinations tend to fall somewhere between the occurrence rates published in Petigura et al.(2013) and 
Burke et al. (2015) for G-dwarfs, which represent a range of about 4 in the warm rocky planet regime, with a tendency to be closer to 
Burke et al. (2015). For example, the geometric mean combination which we use in this paper is about 25% lower than 
Burke et al. (2015) for warm rocky planets, though it is significantly higher than Petigura et al. (2013). It should also be stressed 
that extrapolation is implied when integrating fitted power laws into cold planets or very small planet sizes, so the numbers in 
those regions remain very unreliable.



 To judge the robustness of the SAG13 occurrence rate estimates, independent occurrence rates were calculated 
using the inverse 
detection efficiency method based on the data from the DR25 catalog.
The occurrence rate per bin, η, is given by 

    η   =   1/n_⋆Σ_i^n_p1/ comp_i


Where comp_i is the survey completeness evaluated at the radius and orbital period of each planet in the bin, n_⋆ is the number of stars surveyed, n_p is the number of planets in each bin.
The planet list is taken from <cit.>, using a disposition score cut of 0.9.
The stellar properties are taken from <cit.>, removing giant stars with logg < 4.2
 and culling the G dwarfs by selecting stars with 5300 K < T_eff < 6000 K.
The completeness of each star is calculated with  <cit.>, and the survey completeness is 
calculated by averaging over all stars that were successfully searched for planets according to the  flag.

The above equation assumes that the vetting completeness (the fraction of planet transit signals 
(TCEs, Threshold Crossing Events) properly classified as planet candidates)
 and reliability 
(the fraction of transiting candidates that are not caused by instrumental artifacts or statistical false alarms)
 are 100 %.
The vetting completeness and reliability are very important for small planets, especially at long orbital periods. 
The vetting completeness decreases substantially when one employs a high score cut on the DR25 catalog while the
 reliability approaches 100 %.  The net effect is that occurrence rates are likely to be underestimated by 
ignoring both corrections.
 It should be noted that all SAG13 calculations also ignored vetting completeness and reliability. 
We note that more work needs to be done to 
do a reliable occurrence rate calculation. 


For regions with no planet detections (low instellation, long period orbits), occurrence rates were estimated with a 
parametric function that is a broken power-law in period and radius. Free parameters were constrained using the 
Exoplanet Population Observation Simulator (, Mulders et al. in prep).  generates planet 
populations from this parametrized description using a Monte Carlo simulation, and conducts synthetic observations 
using the survey completeness from the DR25 catalog. The synthetic observable populations are compared with the 
observed planet distribution from Kepler in the range P=[2,400] days and R_p=[0.5,8]  R_⊕ , and the 
posterior parameters are estimated using  (Foreman-Mackey et al. 2013). Binned occurrence rates are 
calculated by marginalizing the posterior parametric distribution over the bin area, and taking the 50% and 16%
 and 84% percentiles for the mean and 1-sigma error, respectively. 

Table <ref> provides the occurrence rates calculated from Eq.(<ref>) for the same bins as in Tables 1 & 2.
The values are more or less consistent within the uncertainties of SAG13 η_basl from Table 2. However, the
extrapolations into the cold planet regimes (low instellation fluxes) results in a disagreement between SAG13
values and from Eq. (<ref>). This is expected, considering that (1) the cold regimes do not have any planet
detections and any extrapolations are expected to wildly deviate (even between the methodologies) from the
true distribution, and (2) the SAG13 rates are a combination of several individual methodoligies. 

[h!]
 Occurrence rates of planets in different boundaries, defined in Table <ref> classification scheme. The η_planet values are estimated using SAG-13 occurrence rates. The ⋆ values are based on extrapolation and 
therefore are very uncertain.



Planet Type (Stellar Flux range)     Radius (R_⊕)     η_pess     η_basl     η_opt
 
 Hot rocky (182-1.0)     0.5-1.0     0.22     0.67     2.04 
 
Warm rocky (1.0-0.28)     0.5-1.0     0.09^⋆     0.30^⋆     1.04^⋆ 
 
Cold rocky (0.28-0.0035)     0.5-1.0     0.50^⋆     1.92^⋆     7.61^⋆
 
Hot super-Earths (187-1.12)     1.0-1.75     0.21     0.47     1.04 
 
Warm super-Earths (1.12-0.30)     1.0-1.75     0.087    0.21     0.54 
 
Cold super-Earths (0.30-0.0030)     1.0-1.75     0.50^⋆     1.42^⋆     4.14^⋆
 

Hot sub-Neptunes (188-1.15)     1.75-3.5     0.29     0.48     0.79
 
Warm sub-Neptunes (1.15-0.32)     1.75-3.5     0.12     0.22     0.41 
 
Cold sub-Neptunes (0.32-0.0.0030)     1.75-3.5     0.77^⋆     1.63^⋆     3.52 ^⋆
 

Hot sub-Jovians (220-1.65)     3.5-6.0     0.05     0.07     0.12 
 
Warm sub-Jovians (1.65-0.45)     3.5-6.0     0.04     0.07     0.13 
 
Cold sub-Jovians (0.45-0.0030)     3.5-6.0     0.58^⋆     1.35^⋆     3.19^⋆ 
 

Hot Jovians (220-1.65)     6.0-14.3     0.028     0.056     0.11 
 
Warm Jovians (1.65-0.40)     6.0-14.3     0.023     0.053     0.12 
 
Cold Jovians (0.40-0.0025)     6.0-14.3     0.34^⋆     1.01^⋆     3.07^⋆ 
 

 





[h!]
 Occurrence rates calculated from Eq.(<ref>). Comparing with the η_basl values 
from Table <ref> from SAG-13 occurrence rates, the SAG13 values are more or less consistent with
the values given below within the uncertainties. We use the η_basl values from Table <ref>
to calculate the exoplanet yield estimates in section 2.2. As with Table <ref>, the ⋆ values
are extrapolated.



Planet Type (Stellar Flux range)     Radius (R_⊕)     η 
 
 Hot rocky (182-1.0)     0.5-1.0     0.552^+0.195_-0.150 
 
      


Warm rocky (1.0-0.28)     0.5-1.0     0.215^+0.148_-0.099^⋆ 
 
      


Cold rocky (0.28-0.0035)     0.5-1.0     1.09^+1.48_-0.755^⋆ 
 
      


Hot super-Earths (187-1.12)     1.0-1.75     0.374^+0.068_-0.056  
 
      


Warm super-Earths (1.12-0.30)     1.0-1.75     0.145^+0.071_-0.061 
 
      


Cold super-Earths (0.30-0.0030)     1.0-1.75     0.78^+0.86_-0.52^⋆ 
 
      


Hot sub-Neptunes (188-1.15)     1.75-3.5     0.356^+0.049_-0.047 
 
      


Warm sub-Neptunes (1.15-0.32)     1.75-3.5     0.147^+0.058_-0.057  
 
      


Cold sub-Neptunes (0.32-0.0.0030)     1.75-3.5     0.85^+0.88_-0.57^⋆  
 
      


Hot sub-Jovians (220-1.65)     3.5-6.0     0.113^+0.019_-0.018  
 
      


Warm sub-Jovians (1.65-0.45)     3.5-6.0     0.051^+0.021_-0.020  
 
      


Cold sub-Jovians (0.45-0.0030)     3.5-6.0     0.279^+0.31_-0.18^⋆  
 
      


Hot Jovians (220-1.65)     6.0-14.3     0.004^+0.011_-0.004  
 
      


Warm Jovians (1.65-0.40)     6.0-14.3     0.002^+0.004_-0.001  
 
      


Cold Jovians (0.40-0.0025)     6.0-14.3     0.008^+0.031_-0.007^⋆  
 
      



 








Integrating the SAG13 parametric model across the bin boundaries defined in Table <ref> gives the occurrence rates in 
Fig <ref>, where we have assumed L = 0.95 L_⊙ and M = 0.965 M_⊙. Each class of planet has an 
occurrence that is a mixture of astrophysical effect and an observational bias. For example, even though it is easier 
to detect giant planets in close orbits, their occurrence rate is comparatively smaller (0.05) than close-in 
sub-Neptune or terrestrial-size planets (∼ 0.48, respectively). The implication is that hot, giant planets are 
likely less in number.  

 The trend in Fig. <ref> indicates that the occurrence rate generally increases, from larger planets to the smaller ones in any particular bin. We can define the inner edge of the habitable zone (HZ) as the boundary where 2o starts 
condensing in a terrestrial planet’s atmosphere, and the outer edge of the habitable zone as 2 condensation 
boundary <cit.>. Within this zone, it appears that the terrestrial size planets have 
higher occurrence rates (0.2 - 0.3) than compared to either Jovians (0.053) or Neptunes (0.07) planets. 
However, it should be noted that the 
occurrence rate of terrestrial planets in this regime is severely restricted by low number statistics. 





 §.§ Mission yield estimates

With the planet categorization scheme and associated occurrence rates described above, we estimated the exoplanet yields for each type of planet using the yield optimization code of <cit.>. Briefly, this code works by simulating the detection of extrasolar planets around nearby stars over the lifetime of a mission.  To do so, it distributes a large number of synthetic planets around each nearby star, sampling all possible orbits and phases consistent with the planet definition, illuminates them with starlight, calculates an exposure time for each planet given a set of assumptions about the instrument and telescope, and determines the fraction of planets that are detectable within a given exposure time (i.e., the “completeness"). The code optimizes the exposure time of each observation, as well as which stars are observed, the number of observations to each star, and the delay time between observations, to maximize the yield for a given type of planet.

We first ran the yield code to define the set of observations that maximized the yield of exoEarth candidates. 
 We adopted the same baseline mission parameters defined in Table 3 of Stark et al. (2015) with exception to the OWA, 
which increased from 15 λ/D to 30 λ/D, the contrast for spectral characterization which improved from 
5× 10^-10 to 1× 10^-10, the spectral resolution which increased from R=50 to R=70, and the SNR required 
for spectral characterization which increased from 5 to 10 per spectral channel.
 We also adopted a new definition for exoEarth candidates.  We distributed exoEarth candidates across the 
<cit.> conservative HZ, which ranges from 0.95 - 1.67 AU for a solar twin. The semi-major axis distribution 
followed the analytic SAG13 occurrence rate fits.  ExoEarth candidates ranged from 0.5 — 1.4 Earth radii, with the lower 
radius limit set by 0.8*(a * (L_⋆)^1/2)^-0.5, a reasonable limit following the work of <cit.>.
  All exoEarth candidates were assigned a flat geometric albedo of 0.2. In this paper, the exoEarths are used solely to 
optimize an observation plan. We do not report on the yield of these exoEarth candidates, focusing instead on the yields 
of the classes of planets defined in section 2 when following such an observation plan.
 We then locked this set of observations in place and re-ran the yield code to 
calculate the yield of each planet type discussed above, simply by changing the planet’s input parameters each time. 
For each planet type, we distributed planet radii and orbital period according to the SAG13 distribution. We assumed 
Lambertian phase functions for all planets.

To calculate the brightness of a given planet, we must also know the planet’s albedo. The actual distribution of exoplanet 
albedos is unknown. So for this study, we simply assigned each planet type a single reasonable albedo. We adopted a 
wavelength-independent geometric albedo of 0.2 for rocky planets and 0.5 for all other planets. 


























To calculate an expected yield, we must also assign each planet type an occurrence rate. Table <ref> lists the occurrence rates obtained from <ref> in each bin of planet size and planet type (hot/warm/cold). The histogram plot in Fig.<ref> visualizes the total scientific impact of the habitable planet candidate survey. The y-axis gives the expected total numbers of exoplanets observed (yields), which are also given by the numbers above the bars. By “expected”, we mean the most probable yield after many trials of an identically executed survey. Three sizes of exoplanets are shown, consistent with Table 1. For each planet size, three incident stellar flux classes are shown: hot (red), warm (dark
blue), and cold (ice blue). The boundaries between the classes correspond to the temperatures where metals, water vapor, and carbon dioxide condense in a planet's atmosphere. The warm bin is not the same as the habitable planet candidate bin,
 as it is likely too generous.
























We also calculated each planet's yield when deviating from the baseline mission parameters. Figs. <ref> & <ref> 
show the sensitivity of each planet's yield to changes in a single mission parameter.  Each yield curve has been 
normalized to unity at the value of the baseline mission.
As expected, the yield of hot planets is more sensitive to the IWA than cold planets, and the yield of cold planets is 
more sensitive to OWA than hot planets.  However, surprisingly, the yield of cold Jupiters is quite sensitive to IWA, 
suggesting that an observation plan optimized for the detection of exoEarths will typically detect cold Jupiters in 
gibbous phase near the IWA.  



We note that in general, larger apertures are less sensitive to changes in 
mission parameters than smaller apertures.









§ DISCUSSION


In this section, we discuss caveats to our classification scheme, and its relevance to future missions that plan to detect and/or characterize extrasolar planets.



 §.§ Caveats to the classification scheme

 The boundaries discussed in earlier sections are made out of the necessity for creating a single classification scheme that applies to all planets, and that can translate current planet obervations into predictions of future planet yields. No single scheme will be able to properly capture the complex interactions between myriad planetary processes in exoplanet atmospheres. An analogy with the habitable zone here is useful: The habitable zone is not a good means of determining the habitability of a single planet. Instead, it is best used to understand how many potentially habitable worlds a given mission may be able to probe for evidence of habitability. Similarly, the classification scheme discussed in this paper is probably not the best means for determining whether a single world has a certain combination of clouds decks or upper atmospheric composition. Ultimately, that should be determined by specific observations of such worlds. What this scheme is useful for is to help understand the diversity, and the number, of worlds that can be expected from future missions and ground-based observatories. 

Even given the above overall caveat, more specific caveats exist. Much of the work upon which this classification scheme is based is relatively new, and represents an area of rapidly evolving thought. This means that the resulting classification scheme will also have to evolve as these theories are better elucidated, tested, and refined. Ultimately, when future exoplanet characterization missions and ground-based observatories become operational, they will provide the tests of the various hypotheses contained in the drawing of the boundaries shown in Fig. 2. This again draws similarities to the habitable zone, which has evolved over the years, and will ultimately be determined by missions. But until such observations are made, this scheme (and the HZ) will enable predictions of how many such worlds we will be able to eventually epxlore in detail.

An example of this rapid evolution is the work by myriad groups on the mass-radius relationship of planets <cit.> Clearly, this area of research is rapidly evolving - and further data on 
this relationship is anticipated with continues ground-based observations of planets around M dwarfs combines with TESS 
observations of the same target. Thus, we anticipate changes to the boundaries selected in this study as these new data and models are incorporated into this line of research. We also anticipate that these data will allow for the study of how these dividing lines are a function of instellation received by the planet. If that is the case, the "horizontal" lines in our classification scheme (Figure 2) would instead be diagonal. We used horizontal lines here as a first-order determination of this scheme and because deriving the slope of those lines empirically or theoretically is well beyond the scope of this work. This is something that may be true more generally - the lines drawn in our scheme may all be diagonal, as we think instellation could impact planet size categories and that planet size could impact the instellations at which various cloud decks form.

A second category of caveat is the lack of consideration of many planetary processes and planet/star properties in this study. Kinetic chemistry, atmospheric circulation, internal heat generation, different star types, and differing bulk planet composition are not explicitly considered in this study. Each of these processes could affect the boundaries considered here.

Our boundaries have been determined at equilibrium. However, the distribution of chemical constituents in planetary 
atmospheres can be strongly affected by the so-called chemical quenching, in which as the material moves, the temperature 
and pressure within the gas mixture change and chemical reactions may become slower, potentially reaching a certain state 
for which the chemical abundances are `frozen' <cit.>. Also, photochemistry could also impact the ability of a 
condensate to form.
This is most evident in the potential for secondary aerosols that form when photochemical byproducts cause supersaturation and condensation in an atmosphere. This process occurs on many solar system worlds <cit.>,
 is thought to have occurred on Archean Earth <cit.>, and likely occurs on exoplanets <cit.>.
 However, these processes themselves  will likely represent a significant overprint over the transitions that are proposed here, where within one of the planet categories we propose, planets that are closer to their host star will have a gradual and increasingly important contribution from non-equilibrium chemistry. 

But others may also exist, including the potential for other photochemical byproducts to form in an atmosphere and be evidenced in the planetary spectrum. However, we expect these will be secondary in importance compared to the optical depth and cold-trapping effects of discrete cloud layers.

The specific properties of the planet and star are also not considered here. Atmospheric circulation, internal heat generation, and the stellar energy distriburion (SED) are also ignored here. These could all impact the specific position of the boundaries we propose. Changing the SED can lead to changes in the planetary albedo, which would affect the position of these boundaries just as they affect the habitable zone <cit.>. Greater internal heat generation, or decreased convection could lead to a suppression of cloud formation. These are not independent variables, as increased heat generation should lead to increased convection and the SED could also impact circulation as it impacts the altitude-dependent deposition of energy in the atmosphere. Although these properties and processes are not included here, they could be in future papers, just as the consideration of the habitable zone now includes consideration of internal heat sources <cit.>, atmospheric circulation and convection. Finally, we do not consider the impact of the bulk chemical composition of the planet. Clearly, a planet that is extremely carbon poor will be less likely to condense  CO_2 or CH_4.

Lastly, we note that the simulations of the scheme in this paper considers a limited parameter space, due to its focus on a near-term problem faced by the exoplanet community. Given the ongoing studies of  flagship missions that propose to directly image and spectroscopically characterize extrasolar planets - HabEx and LUVOIR - we require an ability to discuss and compare the exoplanet yields from such missions. This will also help us understand how these missions will complement the discoveries that will be made by already-planned ground-based and space-based observatories. The focused utility shown in this manuscript comes at the expense of other applications of the proposed classification scheme. It could also be used to simulate the yields from future transit spectroscopy missions such as JWST or CHEOPS. It also helps us understand how detection missions such as TESS, PLATO, and  WFIRST will provide complementary discoveries to past exoplanet detections, and increase the total expected diversity of known planets. These are all worthy applications of the scheme we propose, but we save these applications for future manuscripts.



 §.§ Application of classification scheme to future space-based direct imaging missions


These estimates of the abundances of different planet types allow for projections of the yields of future exoplanet missions. As mentioned in the introduction, two such missions are under study in advance of the next astrophysics Decadal Survey:  HabEx and LUVOIR. Most of the discussion of these missions - and past direct imaging studies - has focused on their ability to find and characterize rocky planets in the habitable zones of nearby stars. However, the observations enabled by such missions would also bring an ability to observe other kinds of worlds. All planets between the inner working angle and outer working angle of the starlight suppression at the time of observation can be observed. And many of these worlds will be brighter than rocky planets in the habitable zone. Therefore, even observation strategies that are optimized for maximizing the yield of rocky planets in the habitable zone will also yield observations of a considerable diversity of worlds. 

None of these yield estimates should be taken as simulations of the yields of the architectures the HabEx and LUVOIR teams are studying. For one, these simulations hold everything except telescope diameter constant, as a way to demonstrate how planet diversity scales with telescope size and as a way to show end-member yields across the range of mission sizes being explored (For a more extensive options of varying quantities other than the telescope diameter, see here: <http://jt-astro.science:5100/multiplanet_vis>). More importantly, the two studies will explore four different architectures which may have different instrument properties, assume different levels of technological development, and different starlight suppression techniques. None of those nuances are considered here. We refer the reader to the eventual reports from these teams for better estimates of the yields from these missions. What we present here should instead be considered an example of how these planet categories can be useful to such yield simulations.

Our simulations predict that a 4m-class mission 


would observe a great diversity of worlds (Figure <ref>). At that scale, a mission whose observations 
are designed to maximize the yield of potential Earths will also yield the detection and characterization of all of the 
planet types discussed here, with the exception of hot Jupiters. Hot Jupiters are not observed by a 4m-class mission 
because the tight inner working angle of the 4m mission, and because of the low abundances of hot Jupiters (see Table 2)
. But a total of up to 90 planets are observed, including up to 12 rocky/super-Earths planets in the habitable zone. 
The diversity of these simulations will allow tests - albeit on a small sample size within each bin - of the planet 
bins proposed here.

An 8m-class mission yields even more planets (Figure <ref>). We predict that class of mission would observe over 
300 planets, including at least a few planets in each of the planet classes proposed in this manuscript.  Although 
neither LUVOIR nor HabEx is explicitly considering an 8m mission, this size sits at the dividing line between the two 
studies and represents an approximation of what a large HabEx or a small LUVOIR could enable. 

The larger versions of LUVOIR (Fig. <ref>) bring the ability to not only observe planets, but to test the occurrence of 
different features 
within each of the planet bins. It would observe dozens of each planet type, providing larger sample sizes which enables to study each 
planet type as a population. As shown by Stark et al. 2014, a sample of 30 planets of a given type would allow us to be sensitive to any
 feature that has at least a 10% chance of occurring. For this class of mission, there would be a total of 
up to 1000 worlds. This includes several hundred
Neptune-size planets with orbits that likely cause water clouds - but not carbon dioxide clouds - to form. The least populated bins are 
for hot Jupiters, which only contains ∼ 15 detections, and for “cold rocky” which would have 2 clouds, 
which contains only few 
detections. But all other bins have at least 50 detections. As mentioned above, 30 planets of a given type would allow 
us to be sensitive 
to any feature that has a ∼ 10% chance of occurrence within that planet type (Stark et al. 2014).
 It is this large number of observation that would allow the myriad hypotheses contained in this manuscript to be 
tested. For example, the presence/absence of various cloud types could be plotted as a function of the energy incident 
at the top of exoplanet atmospheres. And the absence of absorption features associated with cold-trapping could be 
measured in each of these bins. This logic has been applied to habitable planets before <cit.>; here we 
demonstrate that it also applies to planets beyond the habitable zone.











§ CONCLUSIONS



NASA concept mission studies that are currently underway,  like LUVOIR and HabEX, are expected to discover a large 
diversity of exoplanets. These larger missions would provide large enough sample sizes that we could study each planet 
type in the context of its “relatives” in the Solar system or currently known exoplanets.
We present here a classification scheme for exoplanets in planetary radius and stellar flux bins, based on chemical species' condensation sequences in planetary atmospheres. This chemical behavior of gases is dependent on pressure, temperature and metllicity, and we primarily focus on condensation curves of sphalerite ZnS, 2o, 2 and 4. The order of condensation of these species represent the order in which the boundaries of our classification scheme are defined. We then calculated the occurrence rates of different classes of exoplanets within these boundaries, and estimated the direct imaging mission yields for various telescope sizes. While the main focus of future flagship exoplanet direct imaging missions is to characterize a habitable world for bio-signatures, the missions will also have the ability to observe other kinds of planets within the system. Therefore, distinguishing features that separate planet categories based on current observables (planet radius and incident stellar flux), and a scheme to illustrate these categories, is essential in calculating the expected direct imaging mission yields, and correspondingly choosing an optimal observational strategy.


        

The authors would like to thank an anonymous reviewer for their constructive comments. We also thank Aki Roberge,
Mark Marley and Eric Lopez for providing comments on the work and the manuscript. The authors would like to 
gratefully acknowldege the following colleagues who have contributed to SAG-13 community occurrence rate project: 
Chris Burke, Joe Catanzarite, Courtney Dressing, Will Farr, B. J. Fulton, Daniel Foreman-Mackey, Danley Hsu, 
Erik Petigura, Wes Traub.

We would also like to thank the NASA's Kepler 
mission for organizing the “hack week” from October 13 - 15 2015, that provided a forum for project scientists and community 
researchers to interact, and led to significant input into this manuscript.
     R. K and S. D. G. gratefully acknowledge funding from NASA Astrobiology Institute's  Virtual
       Planetary Laboratory lead team, supported by NASA under cooperative agreement
       NNH05ZDA001C. E.H. was supported by an appointment to the NASA Postdoctoral Program at NASA Goddard Space Flight Center, administered by Universities Space Research Association through a contract with NASA.




[Abe et al.(2011)]Abe2011
        Abe, Y., Abe-Ouchi, A., Sleep, N. H. et al. 2011. Astrobiology, 11, 443
        
        [Anglada-Escudé et al.(2016)]Anglada-Escude2016
Anglada-Escudé, G., Amado, P. J., Barnes, J. et al. 2016, Nature, 536, 7617, 437

[Arney et al.(2016)]Arney2016
Arney, G., Domagal-Goldman, S. D., Meadows, V. S. et al. 2016. Astrobiology, 16, 11

[Barnes et al.(2013)]Barnes2013
Barnes, R., Mullins, K., Goldblatt, C. et al. 2013. Astrobiology, 13, 225

[Bonfils et al.(2013)]Bonfils2013
Bonfils, X., Delfosse, X., Udry, S., et al. 2013. , 549, A109

[Burke et al.(2015)]Burke2015
Burke, C. J., et al. 2015. , 809, 8

[Burke & Catanzarite(2017)]BurkeCat2017
Burke, C. J., & Catanzarite, J. Kepler Science document, KSCI-19111-002, Edited by Michael R. Haas and Natalie M. Batalha

[Burrows & Sharp(1999)]BS1999
Burrows, A., & Sharp, C. M. 1999, , 512, 2

[Burrows et al.(2004)]Burrows2004
Burrows, A., Sudarsky, D., & Hubeny, I. 2004. , 609, 1

[Burrows(2005)]Burrows2005
Burrows, A. 2005. Nature, 433, 261

[Cahoy et al.(2010)]Cahoy2010
Cahoy, K. L., Marley, M. S., & Fortney, J. J. 2010, , 724, 1

[Cassan et al.(2012)]Cassan2012
Cassan, A. et al. 2012. Nature, 481, 7380

[Catanzarite & Shao(2011)]CS2011
Catanzarite, J., & Shao, M. 2011. , 738, 151

[Charnay et al.(2015)]Charnay2015
Charnay, B., et al. 2015, , 813, L1

[Chen & Kipping(2017)]ChenKipping2017
Chen, J., & Kipping, D. 2016. , 834, 17

[Coughlin et al.(2016)]Coughlin2016
Coughlin, J. L. et al. 2016. , 224, 12

[Courtin(2005)]Courtin2005
Courtin, R. 2005. Space Science Reviews, 116, 1-2

[Dittman et al.(2017)]Dittman2017
Dittman, J. A., Irwin, J. M., Charbonneau, D. et al. 2017. Nature, 544, 333 

[Dong & Zhu(2013)]DZ2013
Dong, S., & Zhu, Z. 2013. , 778, 53

[Dressing & Charbonneau(2013)]DC2013
Dressing, C., & Charbonneau, D. 2013. , 767, 95

[Dressing & Charbonneau(2015)]DC2015
         Dressing, C., & Charbonneau, D. 2015. , 807, 1
[Evans & Hubbard(1972)]EH1972
Evans, D., & Hubbard, W. 1972. Nature, 240, 162

[Foreman-Mackey et al.(2014)]FM2014
Foreman-Mackey, D., Hogg, D. W., Morton, T. D. 2014. , 795, 1

[Fortney(2005)]Fortney2005
Fortney, J. J. 2005. MNRAS, 364, 469

[Fortney et al.(2007)]Fortney2007
Fortney, J. J., Marley, M. S., Barnes, J. W. 2007. , 659, 2

[Fressin et al.(2013)]Fressin2013
Fressin, F., Torres, G., Charbonneau, D., et al. 2013. , 766, 81

[Fulton et al.(2017)]Fulton2017
Fulton, B., J., Petigura, E. A., Howard, A. W., et al. 2017, , 154, 3

[Gaidos(2013)]Gaidos2013
Gaidos, E. 2013, , 770, 90

[Gillon et al.(2017)]Gillon2017
Gillon, M., Triaud, A. H. M. J., Demory, B. O. et al. 2017. Nature, 542, 456

[Gladstone et al.(2016)]Gladstone2016
Gladstone, G. R., Stern, S. A., Ennico, K. et al. 2016. Science, 351, 6279

[Godolt et al.(2015)]Godolt2015
Godolt, M., Grenfell, J. L., Hamann-Reinus, A. et al. 2015. Planetary and Space Science, 111, 62

[Haqq-Misra & Kopparapu(2015)]HM2015
Haqq-Misra, J., & Kopparapu, R. K. 2015. MNRAS, 446, 428

[Haqq-Misra et al.(2016)]HM2016
Haqq-Misra, J., Kopparapu, R. K., Batalha, N. 2016. , 827, 2

[Howard et al.(2012)]Howard2012
Howard, A. W. et al. 2012, , 201, 15

[Kane et al.(2014)]Kane2014
Kane, S. R., Kopparapu, R. K., Domagal-Goldman, S. 2014. , 794, L5

[Kane et al.(2016)]Kane2016
Kane, S. R., Hill, M. L., Kasting, J. F. et al. 2016. , 830, 1

[Kasting et al.(1993)]Kasting1993
        Kasting, J., F., Whitmire, D., P., & Reynolds. R. T. 1993, Icarus, 101, 108

[Kopparapu et al.(2013)]Kopp2013
        Kopparapu, R. K., Ramirez, R., Kasting, J. F., Eymet, V., Robinson, T. D., Mahadevan, S.,
        Terrien, R. C., Domagal-Goldman, S. D., Meadows, V., & Deshpande, R. 2013, , 765, 131

[Kopparapu(2013)]KoppOccr2013
Kopparapu, R. K. 2013, , 767, L8

   [Kopparapu et al.(2014)]Kopp2014
    Kopparapu, R. K., Ramirez, R. M., SchottelKotte, J. et al. 2014. , 787, L29

[Kopparapu et al.(2016)]Kopp2016
Kopparapu, R. K., Wolf, E. T., Haqq-Misra, J. et al. 2016. , 819, 1

[Kopparapu et al.(2017)]Kopp2017
Kopparapu, R. K., Wolf, E. T., Arney, G. 2017, , 845, 5

[Leconte et al.(2013a)]Leconte2013a
Leconte, J., Forget, F., Charnay, B., Wordsworth, R., & Pottier, A. 2013a. Nature, 504, 268

[Leconte et al.(2015)]Leconte2015
Leconte, J., Wu, H., Menou, K., & Murray, N. 2015. Science, 347, 632

[Lopez & Fortney(2014)]LF2014
Lopez, E. D., & Fortney, J. J. 2014. , 792, 1

[Madhusudhan et al.(2016)]Madhu2016
Madhusudhan, N., Agúndez, M., Moses, J. I., Hu, Y. 2016. Space Science Reviews, 205, 1 

[Marley et al.(2007)]Marley2007
Marley, M. S., Fortney, J., Seager, S., & Barman, T. 2007. Protostars and Planets V, 733

[Mathur et al.(2017)]Mathur2017
Mathur, S., Huber, D., Batalha, N. M., et al. 2017. , 229, 2

[Morley et al.(2012)]Morley2012
Morley, C. V., Fortney, J. J., Marley, M. S. et al. 2012. , 756, 2

[Morley et al.(2013)]Morley2013
Morley, C. V., Fortney, J. J., Kempton, E. M.-R. et al. 2013. , 775, 33

[Morton & Swift(2014)]MS2014
Morton, T. D., & Swift, J. 2014, 791, 10

[Morton et al.(2016)]Morton2016
Morton, T. D. et al. 2016. , 822, 86

[Moses et al.(2011)]Moses2011
Moses, J. I., Visscher, C., Fortney, J. J. et al. 2011. , 737, 15

[Moses et al.(2014)]Moses2014
Moses, J. I. 2014, Nature, 505, 7481

[Mulders et al.(2015)]Mulders2015
Mulders, G. D., Pascucci, I., Apai, D. 2015. , 798, 2

[Parmentier & Guillot(2014)]Parmentier2014
Parmentier, V., & Guillot, T. 2014. , 562, A133

[Parmentier et al.(2015)]Parmentier2015
Parmentier, V. et al. 2015. , 574, A35

[Pierrehumbert & Gaidos(2011)]PG2011
Pierrehumbert, R. T., & Gaidos, E. 2011. , 734, L13

[Petigura et al.(2013)]Petigura2013
Petigura, E. A., Howard, A. W., Marcy, G. W. 2013, PNAS, 110, 48

[Pollack et al.(1996)]Pollack1996
Pollack, J. et al. 1996. Icarus, 124, 1

[Ramirez & Kaltenegger(2017)]RK2017
Ramirez, R. M., & Kaltenegger, L. 2017. , 837. L4

[Robinson & Catling(2012)]RC2012
Robinson, T. D., & Catling, D. C. 2012. , 757, 104

[Robinson & Catling(2014)]RC2014
Robinson, T. D., & Catling, D. C. 2014. Nature geoscience, 7, 1

[Rogers & Seager(2010a)]Rogers2010a
Rogers, L., & Seager, S. 2010a, , 712, 2

[Rogers & Seager(2010b)]Rogers2010b
Rogers, L., & Seager, S. 2010b, , 716, 2

[Rogers(2015)]Rogers2015
Rogers, L. 2015. , 801, 41

[Rowe et al.(2014)]Rowe2014
Rowe, J. F. et al. 2014. , 716, 45

[Silburt et al.(2015)]Silburt2015
Silburt, A., Gaidos, E., Wu, Y. 2015, , 799, 2

[Segura et al.(2003)]Segura2003
Segura, A. et al. 2003. Astrobiology, 3, 4

[Selsis et al.(2007b)]Selsis2007b
Selsis, F. et al. 2007b. , 476, 137

[Shields et al.(2013)]Shields2013
Shields, A. L., Meadows, V. S., Bitz, C. M., et al. 2013. Astrobiology, 13, 8

[Shields et al.(2014)]Shields2014
Shields, A. L., Bitz, C. M., Meadows, V. S. et al. 2014. , 785, L9

[Stark et al.(2014)]Stark2014
Stark, C. C., Roberge, A., Mandell, A., Robinson, T. D. 2014. , 795, 122

[Stark et al.(2015)]Stark2015
Stark, C. C. et al.. 2015. , 808, 149

[Sudarsky et al.(2003)]Sudarsky2003
Sudarsky, D., Burrows, A., Hubeny, I. 2003. , 588, 1121

[Suzuki et al.(2016)]Suzuki2016
Suzuki, D., Bennett, D. P., Sumi, T. 2016. , 833, 145

[Thompson et al.(2017)]Thompson2017
Thompson, S. E., Coughlin, J. L., Kelsey, H. et al. 2017. submitted to , arxiv:1710.06758

[Traub(2012)]Traub2012
Traub, W. A. 2012. , 745, 20

[Waite et al.(2007)]Waite2007
Waite, J. H., Young, D. T., Cravens, T. E. et al. 2007. Science, 316, 5826

[Wakeford & Singh(2015)]WS2015
Wakefor, H. R., & Singh, D. K. 2015. , 573, A122

[Wakeford et al.(2017)]Wakeford2017
Wakeford, H. R., Visscher, C., Lewis, N. K. et al. 2017. MNRAS, 464, 4247

[Way et al.(2015)]Way2015
Way, M. J., Del Genio, A. D., Kelley, M., Aleinov, I., Clune, T. 2015. arxiv:1511.07283

[Weiss & Marcy(2014)]WM2014
Weiss, L., & Marcy, G. 2014. , 783, L6

[Wolf & Toon(2014)]WT2014
Wolf, E., & Toon, E. 2014. Geophysical Research Letters. 41, doi:10.1002/2013GL058376

[Wolf & Toon(2015)]WT2015
Wolf, E., & Toon, E. 2015. JGRA, 120, 5775

[Wolf et al.(2017)]Wolf2017
Wolf, E. T., Shields, A. L., Kopparapu, R. K., et al. 2017. , 837, 107

[Wolfgang et al.(2016)]Wolfgang2016
Wolfgang, A., Rogers, L., Ford, E. B. 2016. , 825, 19

[Wordsworth et al.(2010)]Wordsworth2010
Wordsworth, R., Forget, F., & Eyment, V. 2010, Icarus, 210, 2, 992

[Yang et al.(2013)]Yang2013
Yang, J., Cowan, N. B., & Abbot, D. S. 2013. , 771, L45

[Yang et al.(2014a)]Yang2014a
Yang, J., Gwenael,  B., Fabrycky, D., & Abbot, D.  2014. , 787, L2

[Zahnle & Catling(2017)]ZC2017
Zahnle, K. J., & Catling, D. C. 2017. , 843, 122

[Zsom et al.(2013)]zsom2013
Zsom, A., Seager, S., de Wit, J. et al. 2013, , 778, 2



