A POPULATION-BASED HABITABLE ZONE PERSPECTIVE

Sampling incoming flux and planet radius—We use two different methods to sample the incoming stellar flux and the planet radius. (1) The incoming stellar flux and planet radius of the exoplanet population are drawn from distributions that are uniform in $\mathrm{log}{F}_{\mathrm{in}}$ and $\mathrm{log}{R}_{{\rm{p}}}.$ These two distributions are best suited to studying the general shape and properties of ${p}_{\mathrm{HZ}}.$ If ${F}_{\mathrm{in}}$ and R_p are uniform in log space, it is possible to study ${p}_{\mathrm{HZ}}$ in regions of the parameter space where exoplanets have not been discovered yet (e.g., exoplanets with ${F}_{\mathrm{in}}\lt 0.1{F}_{\oplus }$). (2) The variables ${F}_{\mathrm{in}}$ and R_p are sampled from observed distributions. The second method is used to study the occurrence rate of potentially habitable planets. We use the Dressing & Charbonneau (2013) data set of transiting planets orbiting M dwarfs, rather than similar data sets for sunlike stars (Petigura et al. 2013; Foreman-Mackey et al. 2014). The former sample is more complete to planets with cooler temperatures than Earth, and it does not rely on interpolation in this temperature range to the same extent as for G dwarfs. This is dually true because both the geometric transit likelihood and the signal-to-noise ratio per transit increase for cool planets orbiting M dwarfs, as opposed to G dwarfs.

The data of Dressing & Charbonneau (2013) is recast using Gaussian kernel density estimation. This was necessary for two reasons. The PDF (or occurrence rate) is given as a function of orbital period and planet radius in Dressing & Charbonneau (2013). However, occurrence rates as a function of incoming stellar flux and planet radius are more useful for climate modeling. The PDF values were described on a coarse grid, which is not well suited for random sampling. For these reasons, a Gaussian kernel density estimator is used to characterize the smooth probability density as a function of the incoming stellar flux and planet radius. The joint PDF is illustrated in Figure 2.

Zoom In Zoom Out Reset image size

The PDF of planet mass—The planet radius is an observationally constrained parameter, but the planet mass is unknown for most small exoplanets at habitable distances from the star. Therefore, planet mass PDFs are used to sample the planet mass based on the planet radius. Some planet mass PDFs are illustrated in Figure 3. All PDFs follow a Gaussian distribution around some mean with a standard deviation that is a fraction of the mean. Small rocky planets are distinguished from large planets that can be super-Earths or mini-Neptunes depending on their bulk density. The PDFs can be optimistic, nominal, and pessimistic. The PDFs are optimistic or pessimistic with respect to the occurrence rate of water worlds. The properties of the PDFs are described below. To guide how the PDFs are created, the mass of pure iron, pure silicate, and pure water planets are illustrated as a function of radius R_p (${M}_{i}({R}_{{\rm{p}}}),$ ${M}_{{\rm{s}}}({R}_{{\rm{p}}}),$ ${M}_{w}({R}_{{\rm{p}}}),$ respectively) following Seager et al. (2007).

Zoom In Zoom Out Reset image size

In the optimistic distribution, the planet mass follows a Gaussian distribution if the planet radius is below 1.79 ${R}_{\oplus }.$ The mean of the distribution is ${\bar{M}}_{{\rm{p}}}=({M}_{{\rm{s}}}({R}_{{\rm{p}}})+{M}_{w}({R}_{{\rm{p}}}))/2,$ and the standard deviation can be $\sigma ={\bar{M}}_{{\rm{p}}}/6$ or 3. Here, M_p is not allowed to be smaller than ${M}_{w}({R}_{{\rm{p}}})$ or larger than ${M}_{i}({R}_{{\rm{p}}}).$ Effectively, such a description means that all planets below 1.79 ${R}_{\oplus }$ have a surface (either solid or liquid), and the occurrence rate of pure water planets is significant. If the planet radius is larger than 1.79 ${R}_{\oplus },$ the mean of the distribution is prescribed by an empirical mass–radius relationship (Equation (3) in Weiss & Marcy 2014): ${\bar{M}}_{{\rm{p}}}=2.69{R}_{{\rm{p}}}^{0.93},$ where M_p and R_p are expressed in Earth mass and radius units, respectively. The standard deviation is again $\sigma ={\bar{M}}_{{\rm{p}}}/6$ or 3. Planets are not allowed to be denser than iron. However, the lower mass boundary is given by ${M}_{w}(1.79{R}_{\oplus }).$ The critical density is given by the density of a pure water planet, and all planets denser than the critical density are considered rocky.

It might appear that the optimistic distribution contradicts solar system observations because the inner planets are rocky. However, some moons of Saturn (e.g., Mimas and Thetys) are primarily made of water ice and only small amounts of rock (the density of Thetys is less than that of water). Although these are small moons, they indicate that water-rich worlds could form at larger distances from the star. If such planets have hydrogen-dominated atmospheres, they could harbor liquid water on their surfaces.

The nominal and pessimistic mass distributions differ in three aspects. In the nominal distribution, if the planet radius is below 1.54 ${R}_{\oplus },$ ${\bar{M}}_{{\rm{p}}}={M}_{{\rm{s}}}({R}_{{\rm{p}}}),$ and planets are not allowed to be smaller than $({M}_{w}({R}_{{\rm{p}}})+{M}_{{\rm{s}}}({R}_{{\rm{p}}}))/2.$ In other words, there are no pure water planets, but there exists half water–half silicate rocky planets. The third difference is that the critical density is given by the density of a half water–half silicate planet. In the pessimistic distribution, ${\bar{M}}_{{\rm{p}}}=({M}_{i}({R}_{{\rm{p}}})+{M}_{{\rm{s}}}({R}_{{\rm{p}}}))/2$ below 1.22 ${R}_{\oplus },$ and planets are not allowed to be smaller than ${M}_{{\rm{s}}}({R}_{{\rm{p}}}).$ That is, most rocky planets are assumed to be similar to the rocky planets in our solar system (without a significant water layer). The critical density in this case is given by the density of a pure silicate planet.

Recently, Rogers (2015) showed that most 1.6 ${R}_{\oplus }$ planets are not expected to be rocky (silicate–iron composition), which at first glance would indicate that the optimistic distribution (with all planets below 1.8 ${R}_{\oplus }$ considered to have a surface) is inconsistent with observations. It is important to note that planets drawn from the optimistic distribution are expected to be volatile rich (not silicate and iron rich), so there is no contradiction here. Furthermore, all exoplanets used in the work of Rogers (2015) receive more than 1.1 ${F}_{\oplus }.$ It is likely that more distant exoplanets (with ${F}_{\oplus }\lt 1$) have different properties because such planets might have formed outside the snow line and thus are more volatile rich (similar to the large moons of Jupiter and Saturn). For these reasons, we believe that the optimistic distribution is not inconsistent with observations.

Atmosphere type—If the planet is rocky (its density is larger than the critical density prescribed by the planet mass PDF), the atmosphere type and surface pressure of the planet are sampled. Five different PDFs are considered for the atmosphere type: H₂, N₂, CO₂ dominated atmospheres are equally likely; there are no H₂ dominated atmospheres and N₂ or CO₂ dominated atmospheres are equally likely; there are only N₂ dominated atmospheres; there are only CO₂ dominated atmospheres; or there are only H₂ dominated atmospheres. Some of these PDFs represent extremes. For example, it would be surprising if N₂ or CO₂ dominated atmospheres in the HZ are the exception and H₂ dominated atmospheres are frequent. However, such a PDF cannot be excluded based on our current knowledge.

Surface pressure—Next, a sample is drawn from a surface pressure distribution, and it is important to consider what the range of surface pressures should be. The minimum surface pressure under which water could exist in liquid form is set by the triple-point pressure of water (611 Pa). If the planet's surface pressure is less than 611 Pa, water exists only either in ice or vapor form. Some PDFs are uniform in $\mathrm{log}(P)$ (see Table 1). In this case, the minimum surface pressure is a parameter, and PDFs with ${P}_{\mathrm{min}}$ = 1, 10², 10⁴ Pa are considered. The surface pressure on Mars is 600 Pa for comparison. The minimum surface pressure is set to 0 Pa for lognormal PDFs with modes (most frequently occurring values) of 10⁴, 10⁵, 10⁶, and 10⁷ Pa.

The maximum surface pressure is set by a condition on the total atmospheric mass. The correlation between the atmosphere mass (${m}_{\mathrm{atm}}$) and surface pressure (${P}_{\mathrm{surf}}$) is expressed as

$$\begin{eqnarray}&&{P}_{\mathrm{surf}}=\displaystyle \frac{{m}_{\mathrm{atm}}{g}_{\mathrm{surf}}}{4\pi {R}_{{\rm{p}}}^{2}},\end{eqnarray} \tag{ 1 }$$

where ${g}_{\mathrm{surf}}$ is the surface gravity of the planet. If the atmosphere is N₂ or CO₂ dominated, the maximum surface pressure is such that the atmosphere mass is 1% of the planet's mass. If the atmosphere is H₂ dominated, the maximum surface pressure is set by ${m}_{\mathrm{atm}}={10}^{-5}{M}_{{\rm{p}}}.$

Examples are shown to validate the choices of maximum atmosphere masses and surface pressures. If the atmosphere mass is 1% on an Earth-like planet, the surface pressure is 10⁴ bars. The effective height⁵ of such an atmosphere is only 100 Km higher than the effective height on Earth (assuming a constant 10 Km pressure scale height). Therefore, the relative difference in radius between an Earth-like planet with ${P}_{\mathrm{surf}}$ = 1 bar and ${P}_{\mathrm{surf}}$ = 10⁴ bars is only ∼2%. The typical uncertainty in the measured planet radius is on the order of 10% in the data of Dressing & Charbonneau (2013), mostly due to uncertainty in the stellar radius. Thus, currently we are unable to distinguish an Earth-analog planet from an Earth-like planet with a surface pressure of 10⁴ bars around M dwarfs. The situation is different if Earth hypothetically had an H₂ dominated atmosphere. If ${m}_{\mathrm{atm}}={10}^{-5}{M}_{\oplus },$ the surface pressure is 10 bars. The pressure scale height of H₂ atmospheres is an order of magnitude larger than N₂ atmospheres of the same temperature due to the low mean molecular weight of H₂. If Earth had a 10 bar H₂ atmosphere, the atmosphere's effective height would be ∼200 Km. Although the relative radius difference is larger than in the previous example (∼4%), it is still below the typical uncertainty of Dressing & Charbonneau (2013).

Surface albedo and relative humidity—If the radius of the rocky planet is larger than 0.5 ${R}_{\oplus }$ and the surface pressure is above the critical-point pressure of water, then we proceed to the next step of our habitability investigation. The lower radius limit is set to avoid sampling exoplanets that have low or negative radii: the kernel density estimator can extrapolate to negative-radii exoplanets. Therefore, the lower radius limit is set by the exoplanet with the smallest radius in the data set of Dressing & Charbonneau (2013). In the next step, the surface albedo and relative humidity PDFs are sampled. Three surface albedo PDFs are used: uniform, normal, and lognormal (see Table 1 for more details). The surface albedo is between 0 and 0.5 around M dwarfs, and between 0 and 1 around solar-like stars. The albedo upper limit around M dwarfs is set by the reduced albedo of ice and snow at near-infrared wavelengths where the spectral energy distributions of M dwarfs peak (Joshi & Haberle 2012; Shields et al. 2013).

The relative humidity PDF can be uniform between 1% and 100%, uniform in log space, and lognormal with a mode at 50%. A relative humidity value of 50% is the effective relative humidity in the troposphere of Earth. That is, the global average vertical relative humidity profile of Earth (Manabe & Wetherald 1967) can be approximated with a constant effective relative humidity of 50% to reproduce the global average surface temperature of Earth (Pierrehumbert 2010).

N₂ and CO₂ mixing ratios—If the atmosphere is N₂ or CO₂ dominated, the mixing ratios of these two components are sampled. In an N₂ dominated atmosphere, the mixing ratio of CO₂ is between 0 and 0.5. The PDF of the CO₂ mixing ratio can be uniform, uniform in log space, and lognormal with the mode of $4\times {10}^{-2},$ which is the CO₂ mixing ratio on Earth. Similarly, if the atmosphere is CO₂ dominated, the mixing ratio of N₂ is between 0 and 0.5, and similar PDFs are adapted as before, with the difference that the mode of N₂ in the lognormal distribution is $3.5\times {10}^{-2},$ which is the N₂ mixing ratio in the atmosphere of Venus.

Surface temperature—If the surface pressure is below 10⁷ Pa, we investigate whether the surface temperature is suitable for liquid water. In principle, water can be liquid at pressures above 10⁷ Pa. However, the plausibility of detecting biosignature gases in dense atmospheres decreases with surface pressure (Seager et al. 2013). The reason is that more biomass is needed to produce the same biosignature mixing ratio if the surface pressure is large. If Earth had a surface pressure of 10⁷ Pa, roughly 100 times more biomass would be necessary to reproduce the current atmospheric O₂ and CH₄ mixing ratios, which is unlikely to occur.

The surface temperature is not drawn from a PDF but is calculated from all the parameters sampled up to this point. The reason is that the sampled values unequivocally determine the surface temperature within the limits of a 1D vertical climate model. The incoming stellar flux in combination with the stellar, planetary, and atmosphere parameters allows us to estimate the surface temperature. The method is described in the next section.

Finally, we note that because the surface pressure cutoffs for the atmosphere mass and for the surface temperature check are motivated by observability, ${p}_{\mathrm{HZ}}$ and the HZ occurrence rates should be considered as observable quantities. In principle, water can exist in liquid form at surface pressures higher than 100 bars, but it is potentially very challenging to identify liquid water on the surfaces of such exoplanets.

HUNTER does not perform climate calculations on the fly. Instead, a large set of atmosphere models are precalculated in the relevant range of parameter space and store the results in look-up tables that are interpolated in run time. Such a procedure saves computational time. The look-up tables need to contain the surface temperature as a function of incoming stellar flux and other parameters. The grid dimensions are atmosphere type (three types of atmospheres), surface pressure (between 611 and 10⁷ Pa covered by five grid points), surface gravity (between 2 and 45 m s⁻², six grid points), relative humidity (between 1% and 100%, 10 grid points), surface albedo (five grid points), N₂ and CO₂ mixing ratios (between 0 and 0.5, six grid points), and surface temperature (between 273 K and the boiling point, which is dependent on surface pressure and relative humidity; the grid has a step size of 1 K). The code loops through these grids and generates all-convective temperature–pressure profiles using the climate model of Zsom et al. (2013). The climate model calculates the incoming stellar flux under which the atmospheres are in radiative–convective equilibrium, and the flux values are stored in a look-up table for HUNTER. The assumption of an all-convective profile is valid if one is interested in the surface climate only. This is because the other parts of the atmosphere are largely optically thin and have a small effect on the surface temperature (Pierrehumbert & Gaidos 2011). A modified version of the 1D climate model described in Zsom et al. (2013) is used to generate the look-up tables. It is a globally averaged vertical climate model that includes greenhouse gas absorption, Rayleigh scattering, and collision-induced absorption of N₂, CO₂, and H₂.

HUNTER interpolates the look-up tables to determine the surface temperature at each sampled scenario. If the surface temperature is between 273 K and the boiling point of water, the planet is habitable. Otherwise, the surface temperature is either too large or too small and water is in vapor or ice form, respectively. Because the surface temperature in the tables is restricted to the freezing and boiling points of water, the surface temperature of planets that cannot harbor liquid water is not calculated.

We perform an analysis on atmosphere stability for exoplanets orbiting M dwarfs. We follow the description of Heng & Kopparla (2012) to evaluate whether exoplanet atmospheres in the synthetic population are stable against atmospheric collapse. We assume that the planets are synchronously rotating, and we calculate the radiative and advective timescales of the atmosphere. If the advective timescale is shorter than the radiative timescale, the atmosphere is considered stable. Heng & Kopparla (2012) use the irradiation temperature and the photospheric pressure to calculate the radiative timescale; we use the surface temperature and the surface pressure. The advective timescale is estimated to be ${t}_{\mathrm{conv}}={R}_{{\rm{p}}}/{c}_{{\rm{s}}},$ where c_s is the sound speed at the surface of the exoplanet. In general, we find that atmospheres with low surface pressure and high mean molecular weight are unstable, in agreement with Heng & Kopparla (2012). We do not study the stability of exoplanet atmospheres around Sun-like stars because their rotation periods do not correlate with their orbital periods.

The shape of the average HZ probability as a function of incoming stellar flux and planet radius is illustrated in Figures 4 and 6. The PDFs of planet radius and incoming stellar flux are uniform in log space. Each grid cell contains an average of $\sim {10}^{6}$ synthetic planet scenarios, and only a fraction of those can harbor liquid water on their surfaces. The shape of ${p}_{\mathrm{HZ}}$ is calculated for planets orbiting an M dwarf with an effective temperature of 3500 K (Figure 4, left) and a solar-like star with an effective temperature of 5900 K (Figure 4, right). The other difference between the two simulations is the range of surface albedos. The maximum surface albedo in all PDFs is 0.5 for planets around M dwarfs due to the decreased albedo of snow and ice in the near-infrared. The maximum surface albedo is 1.0 for planets around solar-like stars.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

There are certain regions of parameter space where it is unlikely that planets harbor liquid water because ${p}_{\mathrm{HZ}}$ is zero in all scenarios considered: if a planet is larger than 3 R_⨁ and if a planet receives more than 3 ${F}_{\oplus }.$ In the former case, the planet is expected to be a predominantly gaseous mini-Neptune with no solid surface. In the latter case, the planet is strongly irradiated, and the surface temperature is above the boiling point of water even if the surface albedo is close to the maximum. In other words, the HZ has a sharp inner edge.

Atmospheric stability criteria show that the atmospheres of some exoplanets orbiting M dwarfs are expected to collapse. The average ${p}_{\mathrm{HZ}}$ around M dwarfs with and without atmospheric stability analysis is illustrated in Figure 4. Mostly close-in exoplanets with low surface pressures or high mean molecular weights are expected to collapse. Distant exoplanets that have high surface pressures with a strong greenhouse effect are generally stable against collapse.

Potentially habitable inner-edge planets orbiting M dwarfs are ideal for atmosphere characterization with transmission spectroscopy for several reasons. The ratio of planet to star area is $\sim 1\%$ (Kaltenegger & Traub 2009), so characterization is within reach with the James Webb Space Telescope (JWST). The transit probability is high because the planets are close to the star. Finally, such planets have a short orbital period and thus transit frequently.

The outer edge of the HZ is potentially smooth because one can find a realistic habitable planet scenario even if the incoming stellar flux is 10⁻² ${F}_{\oplus }.$ Atmospheres that have high surface pressures (${P}_{\mathrm{surf}}\gt 10$ bars) and high greenhouse gas concentrations can harbor liquid water on their surfaces if the stellar irradiation is small. The atmospheres of such distant exoplanets can be best characterized by direct imaging in emission due to their distance from the star.

The uncertainty of the HZ probability is large in Figures 4 and 6 (only the average is shown in the figure). For example, the average probability for an Earth-sized planet that receives 1 ${F}_{\oplus }$ is 0.2 and 0.4 around M dwarfs and solar-like stars, respectively. However, the potential range of ${p}_{\mathrm{HZ}}$ is between 0 and 0.7 for these values. The optimistic estimate of 0.7 is uncertain as well due to the limitations of 1D climate modeling (see Section 4.2 for more details).

Zoom In Zoom Out Reset image size

It is interesting to consider the atmospheric properties of planets that are Earth-like (Earth-sized planets that receive 1 ${F}_{\oplus }$) but not habitable. Such planets can be either too hot or too cold to harbor liquid water. If an Earth-like planet has a thick atmosphere dominated by greenhouse gases such as CO₂ or H₂, then the surface temperature is above the boiling point of water. On the other hand, if the atmosphere is dominated by gases that are spectroscopically largely inactive (e.g., N₂) and the atmosphere contains small amounts of greenhouse gases, then the planet is too cold for liquid water (Zsom et al. 2013).

It is notable how similar Figure 4 left and right are. If the HZ limits are expressed in units of incoming stellar flux, variations in the spectral energy distribution of the star and the surface albedo range have a small impact on the surface climate given that the atmospheric properties of exoplanets are currently unknown. In other words, the impact of the stellar energy distribution on the surface climate becomes important only if the exoplanet's atmosphere and geophysical properties are otherwise well known. The slight difference is that the inner edge of the HZ is somewhat closer to solar-like stars than to M dwarfs.

In this section, the occurrence rate of planets bearing liquid water is estimated using the probabilistic HZ description. The PDFs of incoming stellar flux and planet radius are taken from Dressing & Charbonneau (2013), as discussed in Section 2.1. That is, we focus on the occurrence rate of HZ planets around M dwarfs here. We do not show the average HZ probability because it is very similar to Figure 4 left. Small differences are caused because the incoming stellar flux and planet radius phase space are nonuniformly covered in log space (see Figure 2). As a result, some areas of ${p}_{\mathrm{HZ}}$ are not sampled. The marginalized distributions of ${p}_{\mathrm{HZ}}$ with respect to the incoming stellar flux and planet radius are shown in Figure 6. The marginalized distributions allow us to visualize the uncertainty of ${p}_{\mathrm{HZ}}.$ If we want to convert ${p}_{\mathrm{HZ}}$ to occurrence rates, we need to multiply the cumulative ${p}_{\mathrm{HZ}}$ with the total planet occurrence rate of the data set used. The total planet occurrence rate in the data set of Dressing & Charbonneau (2013) is 1.14 planets per star.

The most optimistic estimate on the occurrence rate is 0.3 habitable planets per star. There are more than 4,000 PDF combinations, and thus 4,000 estimates on the HZ occurrence rate. The occurrence rate is maximized if water worlds are frequent (optimistic planet mass PDF). Rocky planet atmospheres tend to be H₂ dominated, and the surface pressure tends to be large, ∼10 bars. Under such conditions, a large fraction of small planets (between 0.5 and 2 ${R}_{\oplus }$) that receive less than 0.6 ${F}_{\oplus }$ could harbor liquid water on the surface. If the incoming stellar flux is larger than 0.6 ${F}_{\oplus },$ the fraction of planets bearing liquid water  gradually declines to zero at around 2 ${F}_{\oplus }.$

The lower limit on the HZ occurrence rate is zero because we do not know how liquid water is partitioned between the atmosphere, surface, and interior on exoplanets. Atmospheric and geophysical processes could potentially remove all water from the surface. Such processes are, e.g., the photodissociation of water molecules and the loss of hydrogen to space, the take-up of liquid water into the mantle, or the incorporation of liquid water into aqueous minerals (e.g., clay). The most pessimistic scenario in the model predicts an occurrence rate of 10⁻⁵ habitable planets per star. In that case, the planet mass follows the pessimistic PDF (ocean planets are rare), the atmosphere is N₂ dominated, and the surface pressure is high, ∼100 bars.

It is up to future observations to constrain the true distributions of planet and atmosphere compositions and surface pressure. Until then, the occurrence rate will remain weakly constrained, and only optimistic upper limits can be derived from the currently available observational data.

The observationally unconstrained random variables are ranked based on their impact on the surface climate and the HZ occurrence rate (see Figure 7). Planetary properties that are ranked high on the list have the strongest impact on the surface climate and habitability, and atmospheric characterization efforts should be dedicated to measuring these properties to best constrain the habitability of individual planets or planet populations.

Zoom In Zoom Out Reset image size

It is possible to rank the random variables because a large number of synthetic populations are available. The first step in ranking is to determine the total average HZ occurrence rate, which is 0.07 planets/star (solid black line in Figure 7). The next step is to determine the average HZ occurrence rate based on those populations that use a given PDF. Because all random variables are represented with several PDFs, there are as many HZ occurrence-rate estimates as PDFs. If the occurrence rates of one random variable strongly scatter around the total occurrence rate, that random variable has a significant impact on the occurrence rate.

The most important planet property that has a strong impact on the surface climate is the surface pressure, followed by the atmosphere type. The surface pressure is important because pressure broadening of absorption lines can greatly increase the greenhouse effect. The surface pressure of rocky planets and moons in the solar system vary greatly: 70 bars on Venus and 0.006 bars on Mars. It is reasonable to assume that the surface pressure on rocky exoplanets is diverse as well, which will have a strong influence on their surface climate. The atmospheric type also has a strong impact on the surface climate. The atmosphere type in the context of the model describes the most abundant component of the atmosphere (i.e., it can be N₂, CO₂, or H₂). Whether the most abundant component is a greenhouse gas or not has strong implications for the surface climate.

We perform a thought experiment in this section to illustrate how additional observational data will help to better constrain the probability that an exoplanet harbors liquid water, and to highlight what will be the greatest challenge in remotely identifying habitable worlds.

Let us assume that we find a transiting planet with a radius of 1 ${R}_{\oplus }$ that receives 1 ${F}_{\oplus }$ insolation orbiting a nearby M dwarf, and the JWST observes its transmission spectrum. Then, we can retrieve the atmospheric composition and the planet mass from the transmission spectrum (de Wit & Seager 2013). We find that the exoplanet has a mass of 1 ${M}_{\oplus },$ a N₂–O₂ dominated atmosphere with 400 ppm CO₂, and some water vapor is also present. That is, the exoplanet appears Earth-like in all aspects. For simplicity, we make the most optimistic assumption that the planet mass and the atmospheric composition are retrieved exactly, without error.

There are three planet properties that could still remain observationally unconstrained: the surface pressure, the relative humidity, and the surface albedo. The surface pressure is unconstrained because we might not be able to distinguish the rocky surface from a cloud deck. Or the surface pressure is larger than a few bars, so the clear atmosphere is optically thick at all wavelengths, and we cannot be sure how deep the atmosphere extends (see Figure 7 of Benneke & Seager 2012). The relative humidity could remain observationally unconstrained for similar reasons. The surface albedo cannot be measured in transmission, only in reflection. And even then, if the surface pressure is too large, light reflected back from the surface might not reach our telescopes.

Under these assumptions, the value of ${p}_{\mathrm{HZ}}$ can be between 10% and 80%. As a reminder, we showed in Section 3.1 that ${p}_{\mathrm{HZ}}$ can be between 0% and 60% for a planet with 1 ${R}_{\oplus }$ that receives 1 ${F}_{\oplus }$ around an M dwarf with an otherwise unknown atmospheric composition and planet mass. In other words, ${p}_{\mathrm{HZ}}$ goes up as expected for an Earth-like atmospheric composition and Earth-like planet mass. However, the uncertainty of ${p}_{\mathrm{HZ}}$ is still large.

The main reason for the large uncertainty is that the surface pressure is observationally unconstrained. To illustrate this, we fix the surface pressure to be 1 bar but leave the relative humidity and the surface albedo unconstrained. Then, the range of ${p}_{\mathrm{HZ}}$ becomes 60–70%. The unconstrained surface albedo and relative humidity still introduce some uncertainty, but the range of ${p}_{\mathrm{HZ}}$ went down significantly.

The thought experiment presented here illustrates that observational and retrieval efforts are paramount to constrain the surface pressure, if we want to identify habitable exoplanets (Misra et al. 2014).

A 1D climate model is used in our description, which is a crude approximation of real atmospheres, but it is computationally efficient and allows us to estimate the surface climate for a broad range of planetary and atmospheric parameters. Recently it has been shown that 1D models overestimate the global average surface temperature compared to 3D global circulation models (GCMs) on synchronously rotating planets with 1 bar surface pressure (Leconte et al. 2013; Yang et al. 2013). The discrepancy is significant: while a 1D model predicts 340 K surface temperature and runaway greenhouse conditions, GCMs calculate a global average surface temperature around and below the freezing point of water. The discrepancy exists because 1D models cannot capture the efficient night-side cooling.

Even though the shortcomings of 1D models are known, it is computationally unfeasible to generate the look-up table with GCMs. As described in the previous section, the table has seven dimensions for each ST considered. If we simulate only three points along each dimension, we need to run 2,187 GCM simulations, which is unprecedented to the best of the author's knowledge. Furthermore, the rotation period (or the order of the spin–orbit resonance for tidally locked planets) and the obliquity would be two additional observationally unconstrained variables for the GCM runs. Thus, almost 20,000 GCM simulations should be performed if only three points along each axis are considered.

We can estimate under what conditions we expect the discrepancy between 1D and 3D models to decrease. The key is to consider when heat redistribution is efficient and thus the day–night temperature contrast is small. It has been shown that the day–night temperature difference on synchronously rotating planets (trapped in a 1:1 spin–orbit resonance) becomes smaller if the surface pressure is larger (Joshi et al. 1997). This is because the atmospheric relaxation timescale linearly increases with surface pressure (Goody & Yung 1989). If the relaxation timescale is large, the atmospheric parcel is transported back to the day side by circulation before it can cool down. The discrepancy is also expected to decrease if planets are trapped in a spin–orbit resonance that is higher than 1:1. A planet with an initial prograde rotation tends to end up in a high-order resonance,⁶ and an initially retrograde planet can reach synchronous rotation (Makarov et al. 2012). A planet trapped in a high-order resonance does not have a permanent night side. Thus the day–night temperature difference is expected to be smaller than on synchronously rotating planets.

As for the impact of 1D modeling on the results, it is expected that the inner edge would be closer to the host star if a GCM or a modified 1D model were used that can treat night-side cooling. Atmospheres with ∼1 bar surface pressures could still be habitable at incoming flux levels where a 1D model predicts too-high surface temperatures.

One caveat of the approach presented here is that I assume the random variables are independent. That is, the value of a random variable does not affect the PDF of other variables. The assumption is justified until our observational knowledge suggests otherwise. It is possible (and straightforward) to include correlations between random variables, but dependencies remain questionable until empirically proven. I highlight two examples here where dependencies might play a role.

It is intuitive to think that the surface pressure should scale with, e.g., the planet radius (Kopparapu et al. 2014) or potentially with the planet's surface area or mass. Solar system objects with a significant atmosphere do not show such scaling, and thus so far we have no evidence to suggest the existence of a trend. If it does exist, it is potentially a weak one with a significant noise. In the future, it might be difficult to find the trend based on observational data. The reason is that the surface pressure of exoplanets can be difficult (but not impossible, Misra et al. 2014) to constrain because the atmosphere might be optically thick at all wavelengths at high surface pressures, or the atmosphere is cloudy and we observe a cloud deck instead of the surface.

The surface albedo is affected by the surface temperature via the ice-albedo feedback. If the surface temperature approaches the freezing point of water, the surface albedo tends to increase because ice and snow are more reflective. The dependency is expected to have a small influence on the results for two reasons: the ice-albedo feedback is weakened around brown dwarfs (Joshi & Haberle 2012; Shields et al. 2013), and the results suggest that the surface albedo PDFs have a weak influence on the HZ occurrence rate (see Figure 7).

The simplest cloud representation was chosen to minimize the number of observationally unconstrained parameters. The radiative effects of clouds are treated as a surface albedo effect. The surface albedo varies between 0 and 0.5 on planets around M dwarfs, in agreement with the ice albedo effect (Joshi & Haberle 2012; Shields et al. 2013), and it is between 0 and 1 on planets around solar-like stars because the surface albedo of a snowball Earth is high around a solar-like star. Clouds are such a complex phenomena that it is likely that a surface albedo approximation remains our best option if exoplanet populations are modeled.

We emphasize that the python implementation of the algorithm (HUNTER) is publicly available. Therefore it can be modified for various cloud treatments in the future.⁷