A Flexible Bayesian Framework for Assessing Habitability with Joint Observational and Model Constraints

The search for a habitable³ Earth-like planet is one of the largest current undertakings in the fields of astronomy and astrobiology. We now have over 3700 confirmed exoplanets, with over 2200 additional planetary candidates. There are several planned space telescope missions such as the Characterising Exoplanets Satellite (CHEOPS; Broeg et al. 2014), the Planetary Transits and Oscillations of Stars satellite (PLATO; Rauer et al. 2014), the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2014), the Wide Field Infrared Survey Telescope (WFIRST; Spergel et al. 2015), the Large Ultraviolet Optical Infrared Surveyor (LUVOIR; Bolcar et al. 2016), the James Webb Space Telescope (JWST; e.g., Barstow & Irwin 2016; Greene et al. 2016), the Origins Space Telescope (OST; Battersby et al. 2018), and the Habitable Exoplanet Observatory (HabEx; Gaudi et al. 2018) that are intended to help discover or further characterize systems (Borucki et al. 2010; Traub 2012; Batalha 2014; Lissauer et al. 2014), and to narrow down the search for habitable worlds. There is also the future possibility of probing the atmospheres of Earth-like exoplanets using ground-based telescopes like the Extremely Large Telescope (ELT; e.g., Hawker & Parry 2019). We now know that exoplanets are relatively abundant in the habitable zones (HZs) of Sun-like and M-type stars (e.g., Tarter et al. 2007; Batalha et al. 2011; Petigura et al. 2013; Dressing & Charbonneau 2015) including small, rocky Earth-like planets in the HZs of their host stars (e.g., Dittmann et al. 2017; Gillon et al. 2017). The question of how we would detect life on another planet has long been considered, including the idea of estimating the likelihood of life existing elsewhere.

Given the sheer number of planetary systems that have been found and are still being discovered, detailed characterization of every system is impractical if not impossible, so it is necessary to prioritize targets for follow up. It makes sense to approach the search for life on another planet from a statistical perspective, which allows us to evaluate the likelihood of detectability progressively as more information is obtained. The most famous example of this type of probabilistic consideration would be the Drake Equation, but other more recent work has also outlined the necessity for a statistical approach for calculating the potential for planetary habitability (e.g., Bean et al. 2017).

The advantage of approaching the complex problem of planetary habitability within a statistical framework is that we can assign probability distributions, regardless of the level at which the important stellar or planetary parameters have been measured. Better observations would mean a narrower and more accurate distribution, but even a poorly characterized system can be broadly constrained by limits imposed by observational data. We can rely on statistical values from large catalogs of measured values for similar systems in order to assess the habitability potential of a particular planetary system of interest for more in-depth, follow-up characterization. For example, Desch et al. (2018) have outlined a specific observational strategy to move toward the prioritization of exoplanet observations. They first focus on observational data that is relatively more easily obtained, and rank planetary systems for more difficult follow-up observations and characterization.

There is a plethora of statistical data available already considered, including exoplanet orbital properties, characteristics, and distributions of planet types/sizes from radial velocity and transit surveys (Kane et al. 2012, 2016; Foreman-Mackey et al. 2014; Kopparapu et al. 2018), and the frequency of planets around stars of specific spectral types (effective temperatures and mass) and/or metallicity values (e.g., Allende Prieto & Lambert 1999; Santos et al. 2003; Ida & Lin 2005; Bean et al. 2006; Adibekyan et al. 2012; Buchhave et al. 2012; Johnson & Li 2012; Bergemann et al. 2014). Additionally, recent work by Bean et al. (2017) discusses a survey of planetary atmospheric compositions for terrestrial planets, for a comparative planetology approach to habitability that utilizes measurable planetary characteristics and a statistical analysis to answer large-scale questions using only a small sample of objects. This method was developed for calculating estimates for the frequency of habitable environments and the existence of life in terms of the exoplanet catalog as a whole, rather than considering planetary systems individually.

In this paper, we consider the inverse problem. It is important to evaluate a planet's habitability history in a probabilistic sense, especially if important characteristics of the system are not known precisely, or at all. From the stellar side of the problem this most prominently includes a host star's mass, composition, and age. Composition may or may not be measured for a particular star. Age is rarely determined to better than a gigayear for low-mass field stars of ages comparable to the Sun, as techniques for determining more precise ages are far more effective for stars younger than 1 Gy (e.g., Soderblom 2010; Silva Aguirre et al. 2015; Almeida-Fernandes & Rocha-Pinto 2018). Given constraints for the properties of a planetary system from observations (if available) or population statistics (if not), and appropriate theoretical models, can we derive a statistical likelihood that a planet will have particular properties? From an astrobiological perspective, we are interested particularly in whether a planet is likely to have detectable life. It is desirable to develop a flexible formalism that can incorporate various constraints as they are contributed by the community.

We implement a Bayesian statistical approach to considering the long-term habitability potential for exoplanetary systems based on joint constraints from modeling and observation. We present an example focused on host star properties. This uses a theoretical model grid of stellar and HZ evolutionary tracks and associated observable distributions for mass and composition of exoplanet host stars and planetary orbital distance. Since Bayesian statistics use mathematical rules of posterior probability to combine observable data and prior information, we produce an inference more precise than gathered from either source alone. It also accommodates disparate sources of information. The utility of this approach is demonstrated by previous work on a Bayesian analysis to understand the probabilities related to emergence of life (Spiegel & Turner 2012), an inherently interdisciplinary enterprise.

High among the problems of considering detectability is the issue of long-term habitability. We need to be able to assess the likelihood that life might not only be present on a particular exoplanet, but also that we would eventually be able to observe it. It is assumed that the ability of life to establish itself and flourish on a planet's surface remains heavily impacted by how long the environment remains stable to even slightly habitable conditions (e.g., McKay 2014; Dong et al. 2018), thus instantaneous habitability is an insufficient consideration. Life must be present for long enough to measurably alter the environment. We choose a time of 2 Gy as the minimum time for a planet to be in the HZ in order for life to possibly be detectable following Truitt et al. (2015) and Truitt & Young (2017) and motivated by, e.g., Brocks et al. (1999), Kopp et al. (2005), and Crowe (2013). This limit can also easily be adjusted without compromising the framework presented here. In this case study we use the host star as a metric for predicting the long-term habitability potential of a given system; this is an important first step, but an incomplete consideration. We recognize that in order to fully assess the habitability potential of a given planetary system, we must consider not only how the characteristics of the parent star will influence the evolution of the HZ, but the nature of the planets (orbital properties, material, atmosphere, etc.) and the conditions of the circumstellar environment (Rushby et al. 2013; Waltham 2017). In this framework, constraints can be relatively easily added to refine predictions as they are developed.

We acknowledge that the relation between an HZ and habitability itself is a complicated concept, and that the conditions for habitability and detectability are difficult to define. Indeed, each might be characterized differently based on whether the associated calculations focus on planetary atmosphere contributions (perhaps related to the search for biosignatures) and the incoming stellar flux, or if approached from a perspective more related to the thermal evolution and the chemical compositions of the planets themselves (Lammer et al. 2009; Güdel et al. 2014). One could consider the effect on habitability depending on whether the planet is a water world or if it is more of an Earth-like planet with continents and plate tectonics (Kite & Ford 2018). The relevant data, though potentially possible to obtain, are challenging to confidently interpret outside of a statistical framework (Seager & Deming 2010; Rein et al. 2014; Catling et al. 2018; Deming & Seager 2017; Schwieterman et al. 2018; Walker et al. 2018), and there are simply a large number of possibly habitable exoplanet candidates.

We are also aware of the likely complications to habitability when we consider the stellar activity of low-mass (M-type) stars, and that planets orbiting in the HZs of these stars are almost necessarily tidally locked (Joshi et al. 1997). These complications are actually drivers of the statistical approach. Even very approximate HZ predictions can place broad constraints on follow-up prioritization, and the HZ model prior can be treated as a module that can be swapped for predictions from an improved model.

In Section 2 we discuss the methodology and the viability of a direct comparison between the data assembled in the Hypatia Catalog and the Kepler database. We also present a Bayesian framework for HZs with four specific cases for different combinations of mass and metallicity priors. In Section 3 we discuss the interpretation of the posterior probability distributions and present several example planetary systems for consideration. In Section 4 we discuss our conclusions and the implications for this approach.

Our overarching goal is to evaluate the relative promise of known planets for hosting detectable life as a way of prioritizing candidates for follow-up observations. The proximate goal of this paper is to address one component by determining the likelihood that a given planet has been in the host star's HZ for at least 2 Gy (denoted CHZ₂), based on measured stellar properties jointly with stellar models. To be clear, CHZ₂ simply means the distances from a star that are in the star's HZ for at least 2 Gy during the star's lifetime, limited to 12 Gy. We use Bayes' theorem with priors derived from observations and outputs from the Tycho code to evaluate the posterior probability that an exoplanet with certain orbital parameters is in the CHZ₂.

The statistical consideration for the distribution of stellar compositions is based on the data assembled in Hypatia (Hinkel & Burger 2017), a catalog comprised of high-resolution stellar abundance measurements. We also implement data from NASA's Exoplanet Archive (Akeson et al. 2013) to understand whether the Kepler exoplanet host stars are well represented by the abundance distribution derived from the Hypatia catalog. Though the measurements included in Hypatia may be inherently biased due to the requirements for obtaining high-resolution abundance measurements, these are stars in the solar neighborhood that give us an idea of the representative distribution of elemental compositions.

We used the stellar evolution code Tycho (Young & Arnett 2005) to create the catalog of evolutionary tracks discussed in detail in Truitt et al. (2015) and Truitt & Young (2017), which are the basis for the model priors discussed here. The database currently contains models between 0.5 and 1.2 solar masses, with metallicities that fall between 0.1 and 1.5 of solar Z-value, so this is the parameter space we analyze here. (Models are also available for non-solar elemental abundance ratios, but we consider only mass and metallicity here.) Tycho outputs information on stellar surface quantities for each time step of the evolution. For this proof of concept, we use those values to define the inner and outer boundaries of the HZ as a function of stellar age, utilizing equations from Kopparapu et al. (2013, 2014), following from Selsis et al. (2007) and Kasting et al. (1993). These prescriptions parameterize the HZ as a function of the host star luminosity and effective temperature, from which we can calculate the associated time-dependent HZ distance for each stellar evolution track. For a given orbital distance from any star we can predict how long and at what stellar age a planet would remain habitable. Thus, for a perfectly characterized star (mass, metallicity, and age known exactly), we can say with certainty whether a planet has been continuously habitable for ≥2 Gy for a given HZ model.

HZ boundaries are very uncertain and model dependent. This is especially true given the use of both 1D and 3D models (e.g., Kopparapu et al. 2013, Leconte et al. 2013). We use the HZ formulation from Kopparapu et al. (2013, 2014) here because they represent combined 1D and 3D model results that present a relatively simple, conservative HZ estimate. For example, the inner HZ may actually be farther in according to some 3D models that require a higher stellar flux for the runaway greenhouse limit (Leconte et al. 2013). Similarly, the outer HZ may be at a larger distance if additional greenhouse gases like H₂, CH₄, and NH₃ are present (Pierrehumbert & Gaidos 2011). Since we are trying to constrain potentially habitable planets, the conservative HZ estimates we use are appropriate since we would rather err on throwing out potentially habitable planets than further polluting our results with planets that are not really in the HZ.

Real habitability is also very complicated and depends on much more than just the star characteristics and the planet's distance from the star, including the planet's rotation, size, mass, composition, and atmosphere. In this work we are trying to use as simple of a framework as reasonably possible to show how a statistical, Bayesian analysis can be made regarding potentially habitable planets. We emphasize that more sophisticated HZ prescriptions can be substituted in the framework presented here, but the main goal of this work is as a proof of concept.

The models are created for fixed stellar masses and compositions, and observations must be used to associate a predicted track with a given star. The observed values will not be precise values. There will be greater or lesser uncertainty. This may reflect unavoidable observational errors for a precise measurement technique all the way to a property that is not measured at all for a given star. It is very difficult to ascertain the ages of low-mass field stars (Soderblom 2010), so we must rely on predictive theoretical models and statistical analyses (e.g., Silva Aguirre et al. 2015; Almeida-Fernandes & Rocha-Pinto 2018). For a star with a given metallicity and mass and unknown age, we can directly derive a posterior probability that a given orbital distance is within the CHZ₂, which is what we do in this case.

To obtain the likelihood P(CHZ₂) that a planet has spent 2 Gy in the HZ, we convolve the probability that a star has a specified value of mass and metallicity, based upon observational constraints, with the theoretical probability of a planet being in the CHZ₂ for the range of M and Z represented by the observational constraints in a Bayesian sense. Below we describe this process for the general case and the following priors: mass M with a Gaussian uncertainty σ_M, metallicity Z with a Gaussian uncertainty σ_Z, and unmeasured metallicity with no assumption about the shape of the metallicity distribution P(Z), and unmeasured metallicity with P(Z) based on the metallicity distribution derived from the Hypatia abundance catalog of nearby stars, as well as combinations of these priors.

2.1. Comparing to Observations: Hypatia and Kepler

We started with data from the Hypatia catalog for 3861 stars with iron abundance measurements within the defined limits of our catalog (0.1 ≤ Fe ≤ 1.5 Z_⊙). Hypatia includes measurements for almost any particular element of interest, and it would be important to consider specific elemental abundance ratios (particularly for oxygen); unfortunately, the Kepler host star data set only provides measurements for metallicity as the iron abundance [Fe/H] and our closest method for comparison is to use the Tycho catalog scaled metallicity Z-value range. We were able to extract a total of 382 host stars from the Kepler exoplanet database that have metallicity measurements between 0.1 and 1.5 Z_⊙.

We used bins from 0.1 to 1.5, centered on those values, with widths of 0.1. Since the measured abundance values in Hypatia have low precision, they tend to fall on the values between bins, like 0.15 or 0.25. For each bin division we took all values at that number and put half of them in the lower bin and half of them in the higher bin. We used the IDL routine KSTWO to compute a Kolmogorov–Smirnov (KS) comparison test (Massey 1951), which outputs the corresponding initial test D-value, D_KS. The KS test compares the distributions (cumulatively adds up the fraction of data points that are less than or equal to a value), with D_KS being the greatest vertical distance between the two data set distributions. The comparison of the differences between the Kepler and Hypatia data sets is shown in Figure 1.

Zoom In Zoom Out Reset image size

To test the significance of the D-values, we used 3861 relevant values from Hypatia and 382 values from Kepler; we wanted to know what the D-value would be if it was Hypatia compared with itself (ultimately the test is to see whether two data sets are representative of the same distribution). We performed 10,000 trials, selecting a random subset of 382 Hypatia values, and did the KS test between the random 382-value subset and the full Hypatia data set. This represents our ideal D-values, since we know that the subset comes directly from the full data set. With these 10,000 D-values, we compute the mean and standard deviation and compare to our original Hypatia–Kepler D-values. We calculate σ = (D_KS–$\left\langle {D}_{\mathrm{KS},\mathrm{trial}}\right\rangle$)/${\sigma }_{{D}_{\mathrm{KS},\mathrm{trial}}}$. Thus, we conclude from these tests that Hypatia and Kepler are representative of non-matching distributions at a level of 0.95σ (see Table 1). We then applied the Anderson–Darling (AD) test (Scholz & Stephens 1987), given by

$$\begin{eqnarray}\begin{array}{rcl}{A}^{2} & = & \displaystyle \frac{n+m-1}{{(n+m)}^{2}}\left({{\rm{\Sigma }}}_{j}\displaystyle \frac{{h}_{j}({\left(n+m){F}_{j}-{{nH}}_{j}\right)}^{2}}{{H}_{j}(n+m-{H}_{j})-(n+m){h}_{j}/4}\right.\\ & & \left.+\displaystyle \frac{1}{m}{{\rm{\Sigma }}}_{j}\displaystyle \frac{{h}_{j}{\left((n+m){G}_{j}-{{nH}}_{j}\right)}^{2}}{{H}_{j}(n+m-{H}_{j})-(n+m){h}_{j}/4}\right),\end{array}\end{eqnarray} \tag{ 1 }$$

where j runs over the Fe values (low-to-high), F is the Hypatia set, n is the number of values in F, G is the Kepler set, m is the number of values in G, F_j represents the number of Hypatia values less than j plus one-half of those equal to j, G_j is the same but with Kepler, H_j is the same but with both data sets combined, and h_j is the number of values equal to j from the combined data sets. We compute this value and do 10,000 more trials comparing a Hypatia subset to itself. A² represents non-matching distributions at a level of −0.0479336σ.

Table 1.  Comparing the Distributions of the Hypatia and Kepler Data Sets

	KS-Test (D-values)	AD-Test (A2-values)
Initial values	0.0656673	1.11140
10,000 trials, avg	0.0518119	1.15103
10,000 trials, sig	0.0145691	0.826723
Final distribution	0.951011σ	−0.0479336σ

Download table as:  ASCIITypeset image

For the normal standard deviation definition, these values can be translated to mean that the KS test gives a 34% chance that these two data sets are representative of the same distribution, and the AD-test gives a 96% chance. While these two results differ significantly neither of them can confidently rule out the possibility that the Kepler and Hypatia data sets are representative of the same distribution. We conclude from these tests that the Hypatia catalog is a fairly good representation of the Fe distribution of the Kepler planet hosts, and thus we can move forward in comparing these two data sets with our catalog of stellar models.

2.2. A Bayesian Framework for the Habitable Zone

We use Bayes' Theorem to create a statistical framework for assessing long-term planetary habitability. Ultimately, we can calculate conditional probabilities based on a mathematical formulation of observables, in order to estimate the habitability potential for a given system. We need two mathematical functions in order to define this framework: a function for the likelihood of a certain outcome, and a distribution of defined priors. The product of these will give us the posterior distribution. Bayes' Theorem is commonly written as

$$\begin{eqnarray}&&P(A| B)=\displaystyle \frac{P(B| A)\times P(A)}{P(B)},\end{eqnarray} \tag{ 2 }$$

where $P(A| B)$ is the posterior probability of A occurring given that B is true, $P(B| A)$ is the probability of B occurring given that A is true, P(A) is probability of A, and P(B) is the probability of B. Given two possible outcomes, A₁ and A₂, the denominator of Equation (2) can be rewritten using relation $P(B)=P(B| {A}_{1})\times P({A}_{1})+P(B| {A}_{2})\times P({A}_{2})$. We therefore have

$$\begin{eqnarray}&&P({A}_{1}| B)=\displaystyle \frac{P(B| {A}_{1})\times P({A}_{1})}{P(B| {A}_{1})\times P({A}_{1})+P(B| {A}_{2})\times P({A}_{2})}.\end{eqnarray} \tag{ 3 }$$

We apply Equation (3) to Tycho models and observational distributions of stellar properties to calculate the Bayesian posterior probability that an exoplanet with certain orbital parameters is in the CHZ₂. We limit the stellar models to a total age of 12 Gy to take into account the age of the universe. We perform this analysis for four cases, first using B = Z where Z is the metallicity taken from Hypatia, then using B = Z where Z is a measured stellar metallicity with an assumed Gaussian error distribution, next using B = Z, M where Z is metallicity from Hypatia, and mass M is a measured value with an assumed Gaussian error, and finally B = Z, M, where both Z and M are measured values with Gaussian errors.

2.2.1. Case 1: Using a Generalized Stellar Metallicity Distribution

Here, we use the metallicity Z as the only prior in our Bayes' equation and assume a fixed mass, so from Equation (3) we have A₁ = CHZ₂, ${A}_{2}=\,\mathrm{not}\,{\mathrm{CHZ}}_{2}$, and B = Z. A₁ is the possible outcome where a given radius has been in the CHZ for 2 Gy, and A₂ is the outcome wherein the radius is not in the CHZ₂. Equation (3) therefore becomes

$$\begin{eqnarray}&&\begin{array}{l}P({\mathrm{CHZ}}_{2}| Z)\\ =\,\displaystyle \frac{P(Z| {\mathrm{CHZ}}_{2})P({\mathrm{CHZ}}_{2})}{P(Z| {\mathrm{CHZ}}_{2})P({\mathrm{CHZ}}_{2})+P(Z| \mathrm{not}\,{\mathrm{CHZ}}_{2})P(\mathrm{not}\,{\mathrm{CHZ}}_{2})}.\end{array}\end{eqnarray} \tag{ 4 }$$

To obtain the observational Z probabilities, we use the distribution from the Hypatia catalog of nearby stars. For Z we use the iron abundance [Fe/H] values, where Z/Z_⊙ = 10^[Fe/H]. For each calculation we use a fixed mass M_j where the index j refers to the specific mass model used, and only the Tycho model outputs for that mass, linearly interpolated if necessary. We describe how we compute each component of Equation (4) for a desired metallicity Z_j and the fixed mass M_j, for any given radius from the star. We compute P(CHZ₂), the probability that a given radius around a star of mass M_j is in the CHZ₂, by

$$\begin{eqnarray}&&P({\mathrm{CHZ}}_{2})=\displaystyle \frac{{\displaystyle \sum }_{i}{t}_{{\mathrm{CHZ}}_{2},i}{\rm{\Delta }}{Z}_{i}}{{\displaystyle \sum }_{i}{t}_{\mathrm{tot},i}{\rm{\Delta }}{Z}_{i}}\end{eqnarray} \tag{ 5 }$$

where the index i runs through the Tycho model metallicities from Z = 0.1–1.5 Z_⊙ , ${t}_{{\mathrm{CHZ}}_{2},i}$ is the total time the radius is in the CHZ₂ for a given (Z_i, M_j) Tycho model, ${t}_{\mathrm{tot},i}$ is total lifetime (with a maximum of 12 Gy) for the given (Z_i, M_j) Tycho model, and ΔZ_i is Z_i+1 − Z_i, the distance between Z values used. The value of $P(\mathrm{not}\,{\mathrm{CHZ}}_{2})$ is given by

$$\begin{eqnarray}&&P(\mathrm{not}\,{\mathrm{CHZ}}_{2})=1-P({\mathrm{CHZ}}_{2})\end{eqnarray} \tag{ 6 }$$

and the probability that a star has metallicity Z_j if a given radius is in the CHZ₂ is

$$\begin{eqnarray}&&P({Z}_{j}| {\mathrm{CHZ}}_{2})=\displaystyle \frac{{t}_{{\mathrm{CHZ}}_{2},j}\times P({Z}_{j})}{{\displaystyle \sum }_{i}({t}_{{\mathrm{CHZ}}_{2},i}\times P({Z}_{i}))}.\end{eqnarray} \tag{ 7 }$$

P(Z) is a Gaussian fit to the fraction of stars in Hypatia with metallicity Z, given by

$$\begin{eqnarray}&&P({Z}_{i})=0.0937309\times {e}^{-{\left(({Z}_{i}-0.927451)/0.498708\right)}^{2}/2}.\end{eqnarray} \tag{ 8 }$$

We similarly compute $P({Z}_{j}| \mathrm{not}\,{\mathrm{CHZ}}_{2})$, the probability that a star has metallicity Z_j if a given radius around the star is not in the CHZ₂, by

$$\begin{eqnarray}&&P({Z}_{j}| \mathrm{not}\,{\mathrm{CHZ}}_{2})=\displaystyle \frac{({t}_{\mathrm{tot},j}-{t}_{{\mathrm{CHZ}}_{2},j})\times P({Z}_{j})}{{\displaystyle \sum }_{i}(({t}_{\mathrm{tot},i}-{t}_{{\mathrm{CHZ}}_{2},i})\times P({Z}_{i}))}.\end{eqnarray} \tag{ 9 }$$

When calculating the Bayesian posterior probability from Equation (4), we combine Equations (5)–(7), and (9) to find $P({\mathrm{CHZ}}_{2}| {Z}_{j})$. This gives, with a factor of P(Z_j) canceling out,

$$\begin{eqnarray}&&\begin{array}{l}P({\mathrm{CHZ}}_{2}| {Z}_{j})\\ =\,\displaystyle \frac{\tfrac{{t}_{{\mathrm{CHZ}}_{2},j}}{{\displaystyle \sum }_{i}({t}_{{\mathrm{CHZ}}_{2},i}P({Z}_{i}))}\tfrac{{\displaystyle \sum }_{i}{t}_{{\mathrm{CHZ}}_{2},i}}{{\displaystyle \sum }_{i}{t}_{\mathrm{tot},i}}}{\tfrac{{t}_{{\mathrm{CHZ}}_{2},j}}{{\displaystyle \sum }_{i}({t}_{{\mathrm{CHZ}}_{2},i}P({Z}_{i}))}\tfrac{{\displaystyle \sum }_{i}{t}_{{\mathrm{CHZ}}_{2},i}}{{\displaystyle \sum }_{i}{t}_{\mathrm{tot},i}}+\tfrac{{t}_{\mathrm{tot},j}-{t}_{{\mathrm{CHZ}}_{2},j}}{{\displaystyle \sum }_{i}(({t}_{\mathrm{tot},i}-{t}_{{\mathrm{CHZ}}_{2},i})P({Z}_{i}))}\left(1-\tfrac{{\displaystyle \sum }_{i}{t}_{{\mathrm{CHZ}}_{2},i}}{{\displaystyle \sum }_{i}{t}_{\mathrm{tot},i}}\right)}.\end{array}\end{eqnarray} \tag{ 10 }$$

We now have a Bayesian posterior probability at each radius for a given Z and M. We can marginalize over Z in $P({\mathrm{CHZ}}_{2}| Z)$ for a given mass and make a contour plot of all masses, shown in Figures 2–4, and additionally overplot exoplanets with known stellar masses and orbital distances from the Kepler exoplanet archive, shown in the right panel of Figure 4. As an exercise, we can also apply this process to a few known planets, with the results summarized in Table 2 and shown in Figures 5 and 6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2.  Results of Cases 1–4 for Several Solar System Planets and Exoplanets of Interest

Planet	Distance	Mass	σM	Metallicity	σZ	Radius	$\left\langle P({\mathrm{CHZ}}_{2}| Z)\right\rangle$	$\left\langle P({\mathrm{CHZ}}_{2}| Z)\right\rangle$	$\left\langle P({\mathrm{CHZ}}_{2}| Z,M)\right\rangle$	$\left\langle P({\mathrm{CHZ}}_{2}| Z,M)\right\rangle$
	(ly)	(M⊙)	(+/−)	(Z/Z⊙)	(+/−)	(au)	(Hypatia-Z)	(measure-Z)	(Hypatia-Z)	(measure-Z)
							Case 1	Case 2	Case 3	Case 4
Venus	⋯	1.0	0.001	1.0	0.001	0.72	0.000	0.000	0.000	0.000
Earth	⋯	1.0	0.001	1.0	0.001	1.00	0.355	0.322	0.380	0.322
Mars	⋯	1.0	0.001	1.0	0.001	1.52	0.802	0.791	0.810	0.791
Tau Ceti e	12	0.78	0.012	0.28	0.05	0.52–0.58	0.000	0.000	0.000	0.000
Tau Ceti f	12	0.78	0.012	0.28	0.05	1.29–1.35	0.064	0.054	0.194	0.370
Kepler-186f	561	0.54	0.02	0.55	0.12	0.42–0.44	0.837	0.835	0.831	0.835
Kepler-62f	1200	0.69	0.02	0.43	0.04	0.71–0.73	0.529	0.787	0.780	0.834

Note. The radius ranges indicated here are representative of ±1σ semimajor axis orbits, and $\left\langle P({\mathrm{CHZ}}_{2}| X)\right\rangle$ is a marginalization of $P({\mathrm{CHZ}}_{2}| X)$ over the entire radius span.

Download table as:  ASCIITypeset image

2.2.2. Case 2: We Have a Stellar Metallicity Measurement

If we know the Z measurement in the form of Gaussian errors, as the Z' ± σ_Z generally quoted in abundance surveys, we can use this as a more accurate representation of the term P(Z) by using the Gaussian probability

$$\begin{eqnarray}&&P(Z)=\displaystyle \frac{1}{\sqrt{2\pi }}{e}^{\tfrac{1}{2}\tfrac{{\rm{\Delta }}^{\prime} {Z}^{2}}{{\sigma }_{Z}^{2}}},\end{eqnarray} \tag{ 11 }$$

where ${\rm{\Delta }}^{\prime} Z$ is defined as Z' − Z. This P(Z) is used in Equation (10). We apply this method to our planets of interest, with results in Table 2 and Figure 7. For the Sun-related calculations that use measurements and corresponding errors for the M and Z values, we have simply assumed very small errors of 0.001 for both M (M/M_⊙) and Z (Z/Z_⊙).

Zoom In Zoom Out Reset image size

2.2.3. Case 3: Mass Measurement but No Metallicity Measurement

If we know the mass measurement in the form of Gaussian errors, like $M^{\prime} \pm {\sigma }_{M}$, we can, like before, use this as a representation of the term P(M) using the Gaussian probability

$$\begin{eqnarray}&&P(M)=\displaystyle \frac{1}{\sqrt{2\pi }}{e}^{\tfrac{1}{2}\tfrac{{\rm{\Delta }}^{\prime} {M}^{2}}{{\sigma }_{M}^{2}}}\end{eqnarray} \tag{ 12 }$$

where Δ'M is defined as M' − M. This P(M) is then used in Equation (13). If we use mass as another known property in our Bayes' equation, and if we also assume that the stellar mass and metallicity values are independent, Equation (10) becomes

$$\begin{eqnarray}&&\begin{array}{l}P({\mathrm{CHZ}}_{2}| {Z}_{j},{M}_{j})\,=\,\displaystyle \frac{\tfrac{{t}_{{\mathrm{CHZ}}_{2},j}}{{\displaystyle \sum }_{i}({t}_{{\mathrm{CHZ}}_{2},i}P({Z}_{i})P({M}_{i}))}\tfrac{{\displaystyle \sum }_{i}{t}_{{\mathrm{CHZ}}_{2},i}}{{\displaystyle \sum }_{i}{t}_{\mathrm{tot},i}}}{\tfrac{{t}_{{\mathrm{CHZ}}_{2},j}}{{\displaystyle \sum }_{i}({t}_{{\mathrm{CHZ}}_{2},i}P({Z}_{i})P({M}_{i}))}\tfrac{{\displaystyle \sum }_{i}{t}_{{\mathrm{CHZ}}_{2},i}}{{\displaystyle \sum }_{i}{t}_{\mathrm{tot},i}}+\tfrac{{t}_{\mathrm{tot},j}-{t}_{{\mathrm{CHZ}}_{2},j}}{{\displaystyle \sum }_{i}(({t}_{\mathrm{tot},i}-{t}_{{\mathrm{CHZ}}_{2},i})P({Z}_{i})P({M}_{i}))}\left(1-\tfrac{{\displaystyle \sum }_{i}{t}_{{\mathrm{CHZ}}_{2},i}}{{\displaystyle \sum }_{i}{t}_{\mathrm{tot},i}}\right)}.\end{array}\end{eqnarray} \tag{ 13 }$$

For P(Z) we use the Gaussian fit of the fraction of stars in Hypatia with the metallicity Z given by Equation (8). We apply this method to the same planets of interest as previously discussed, with the results summarized in Table 2 and shown in Figure 8.

Zoom In Zoom Out Reset image size

2.2.4. Case 4: Mass Measurement and a Metallicity Measurement

If we know both the mass and metallicity measurements for a star, we can use Equations (11) and (12) along with Equation (13). We can also apply this method to the same planets as discussed previously, with results summarized in Table 2 and in Figure 9.

Zoom In Zoom Out Reset image size

To begin, we discuss the Hypatia-weighted model probabilities for P(CHZ₂). This corresponds to Case 1, which assumes we only know the stellar mass. This provides a quick estimate for the habitability potential of any given system without follow-up observations of composition, as is frequently the case with more distant Kepler candidates that do not have high resolution spectroscopy. The left panel of Figure 2 shows the posterior probability of a given orbital distance having been continuously habitable for at least 2 Gy for metallicity in steps of 0.05 Z/Z_⊙ from 0.1 to 1.5 Z/Z_⊙, assuming a fixed mass of M = 1 M_⊙ and an unconstrained age. These distributions are mapped out in an associated contour plot (right panel), where P$({{\rm{CHZ}}}_{2}| Z)$ is indicated by the color bar at the top. The inner and outer HZ limits at the zero-age main sequence (ZAMS) and the terminal-age main sequence (TAMS) are also included, where the inner boundary is defined by the conservative runaway greenhouse case, and the outer boundary is defined by the maximum greenhouse case, discussed in detail in Kopparapu et al. (2013). For the high Z-cases, the HZ distance is closer-in and narrower, while the highest likelihood of a 2 Gy continuously habitable orbit at some point during the MS evolution occurs between 0.8 and 1.2 Z/Z_⊙, where the Hypatia distribution peaks. Based on observational data, it appears that while low-Z stars have larger HZs, the smaller lifetimes and observational likelihoods give higher-Z stars larger probabilities. This trend is consistent for all masses in Figure 3.

Figure 3 shows the P(CHZ₂) contour plots as a function of mass for each Z-value. For each panel in the figure, the distance from the star is shown on the x-axis, the stellar mass is given on the y-axis, and the color contour illustrates the associated calculated P-value as a visual representation of the habitability potential. As expected for the shorter lifetimes of higher mass stars, the regions of high P(CHZ₂) decrease as mass increases, while the total width of the range of habitable orbits increases due to the faster evolution of the stars.

Figure 4 (left panel) shows the $P({\mathrm{CHZ}}_{2}| Z)$ contour plot for all masses in our model grid (0.5–1.2 M_⊙) marginalized over all Z values (0.1–1.5 Z_⊙) referenced in Figure 3. In the right panel, we overplot Kepler exoplanets whose host stars fall within the boundaries of our parameter space. Since the known Kepler exoplanet host star metallicity distribution is consistent with the Hypatia distribution (Figure 1), this is an appropriate prior.

Figure 5 expands on the results from Figures 2–4. It is representative of the data for the specific exoplanets of interest indicated in Table 2, now overplotted on the corresponding mass and metallicity contour plot (as demonstrated in Figure 3), interpolated to the measured Z-value of the host star. The black dots that indicate the location of the planet on this grid indicate the uncertainty in the orbital distance measurements for the planets. Figures 6–9 (Cases 1–4, respectively) utilize example exoplanetary systems to demonstrate the use of the Bayesian statistical framework for assessing the detectability potential of a given system; we show the Bayesian posterior probability distributions of the habitability potential for several well-known exoplanets of interest: Tau Ceti (e) and (f) (upper left; Tuomi et al. 2013), Kepler-186f (upper right; Torres et al. 2015), Kepler-62f (lower left; Borucki et al. 2013), and the Sun, for comparison. The results for the Sun (shown in the bottom right panels of each of these figures) have particularly interesting implications, in that if we are using our own solar system as a standard for looking at other potentially habitable systems, we might choose to focus on different planets based on the measured values available.

Some of the subpanel plots in Figures 6–9 seem to reach a maximum value plateau around a value of 0.833, which corresponds to a fractional value of 10/12. This maximum value is inherently defined by our method, due to our inclusion of the 12 Gy MS lifetime cutoff. Because we have limited the MS lifetime of each stellar model in our grid to 12 Gy total, and we are looking specifically at the viability of a 2 Gy continuously HZ, the longest possible time that any planet could be in the CHZ₂ is 10 Gy. Thus, the highest possible posterior probability is roughly at the value of 10/12. So, a P-value of 0.835 indicates that the orbital radius of that particular planet would be considered to be continuously habitable for the entire range of the MS lifetime considered (up to the maximum value of 12 Gy).

Figure 6 is for Case 1, where we calculate $P({\mathrm{CHZ}}_{2}| Z)$ using Hypatia and a fixed mass; Figure 7 is for Case 2, where we calculate $P({\mathrm{CHZ}}_{2}| Z)$ using Z measurements and a fixed mass; Figure 8 is for Case 3, where we calculate $P({\mathrm{CHZ}}_{2}| Z,M)$ using Hypatia and a mass measurement; and Figure 9 is for Case 4, where we calculate $P({\mathrm{CHZ}}_{2}| Z,M)$ using Z measurements and a mass measurement. It will be extremely important to obtain the best possible measurements for stellar mass and metallicity to calculate the most accurate P-value for habitability potential. For purposes of comparison in Cases 2–4 we assign a Gaussian uncertainty of σ_M and σ_Z = 10⁻³ for both solar mass and metallicity. The data in the upper right corner of each panel are the associated P-values for the planets and exoplanets of interest.

For the M measurement and Hypatia Bayesian distribution (Case 3; Figure 8), we find the posterior probability of Venus to be in the CHZ₂ is P = 0 (this is expected, as Venus would not be considered habitable by most definitions), Earth is P = 0.380, and Mars is the highest with a P = 0.810. This by itself is due to the fact that in these HZ models, Mars would be continuously habitable over a longer fraction of the Sun's overall MS lifetime. When we impose measurements for both Z and M (Case 4; Figure 9) we again calculate a P-value of P = 0 for Venus, P = 0.322 for Earth, and P = 0.791 for Mars. This is a useful illustration of how the posterior probability distribution will change as constraints are added. With a wide range of possible compositions, there are more evolutionary histories that allow a given orbit to be in the CHZ₂. As a result, the distribution is wider, with a higher posterior probability that a given orbit has been habitable. When the value of Z is more constrained, the distribution is narrower.

The Kepler exoplanets included here have a high likelihood to be terrestrial, based on the Earth similarity index (ESI) discussed in Schulze-Makuch et al. (2011). All have a radius between 0.5 and 1.5 R_⊕ and masses between 0.1 and 5 M_⊕; they also fall within the limits of the conservative HZs around the associated host star. $P({\mathrm{CHZ}}_{2}| Z)$ and $P({\mathrm{CHZ}}_{2}| Z,M)$ values are calculated for each exoplanet of interest. The data in the upper right corner of each panel are the stellar mass, the metallicity value (Z/Z_⊙) of the host star, and the P-value marginalized over the entire orbital range.

The Tau Ceti exoplanets were chosen specifically as a follow-up investigation to the work done in Pagano et al. (2015). In Figure 9, $P({\mathrm{CHZ}}_{2}| Z,M)$ for planet (e) of the Tau Ceti system (blue) is P = 0, since the planet is quite close-in to the star. We consider planet (f) as more solidly dwelling in the HZ of Tau Ceti (red) with P = 0.370, but it still has a relatively low posterior probability for long-term habitability (in terms of the CHZ₂) due to the fact that (f) has likely only recently entered the HZ. These calculations confirm the result of Pagano et al. (2015), who noted that planet (e) would be considered to be in the HZ only with very optimistic assumptions. Planet (f) would initially appear to be a somewhat more robust candidate as a potentially habitable exoplanet, but the associated stellar evolution track shows that due to the increased luminosity of Tau Ceti over its MS lifetime, it is likely the planet has only recently moved into the star's conservative HZ. We find that Kepler-186f has the highest posterior probability, with P = 0.835 and would likely be the best candidate for in-depth follow-up observations with JWST or other applicable missions. Kepler-62f has only marginally less habitability potential, with P = 0.834, and would likely also make a good candidate for follow-up characterization. These posterior probability values are consistent with what we would expect for the overall habitability potential for these planets, in terms of the ESI of Kepler exoplanets, as well as the assessment outlined in Pagano et al. (2015).

The intent of this paper has been to demonstrate that we can use our catalog of Tycho stellar evolution tracks, along with observational data, to estimate the likelihood that an observed exoplanet is in its host star's 2 Gy continuously HZ (CHZ₂) at the time of observation, even if certain stellar properties are poorly measured or unknown. We have shown that the Hypatia distribution of iron abundance measurements in the included stars is a reasonable representation for the metallicities cataloged in the Kepler exoplanet host star database. By constructing a Bayesian statistical framework wherein we can utilize model constraints and observations of individual systems or populations as priors, we can draw meaningful conclusions about the potential for long-term planetary habitability. We have shown that we have the capability to do this for specific well-studied exoplanets of interest. The instantaneous present-day location of the HZ is an insufficient criterion for assessing the detectability of life on a planet. We conservatively assess the 2 Gy continuously HZ, which is based on the timescale for the concurrent evolution of life and conceivably detectable gases in Earth's atmosphere.

A statistical approach toward the characterization of detectable, potentially habitable exoplanetary systems is an important aspect to be considered in the ongoing search for life elsewhere. In this proof of concept, we confirm that stellar mass and composition (specifically, metallicity [Fe/H], or Z/Z_⊙ values) play important roles in the evolution of a star's HZ, as well as influencing an orbiting exoplanet's long-term habitability potential. We have shown that by implementing a Bayesian statistical approach, we can use Tycho stellar evolution models to inform what kind of stars and what specific orbital radii (habitable distances) are potentially the best to look at with follow-up observations with upcoming space telescope missions. Since our catalog of Tycho models is a reasonable representation of the intrinsic spread in values observed for the mass and metallicity of real, nearby host star candidates, we are now better able to leverage the database of stellar evolution tracks to assess the long-term habitability potential of a given planetary system of interest.

The results of this work can also be complemented by deeper, more specific analysis and observations to better understand the potential habitability of a given planet. For example, it is understood that a planet's habitability can vary significantly within the HZ (e.g., Van Laerhoven et al. 2014). When looking at a subset of promising targets, additional planetary and stellar characteristics can be taken into account when trying to better understand the potential for habitability, especially for candidates near the extremes of the HZ.

This kind of statistical assessment gives us more analytical power, since we do not need to know every single stellar or planetary parameter in order to calculate a Bayesian posterior probability and ultimately gauge the habitability potential of a particular stellar system. Conversely, it is relatively straightforward to add additional constraints that increase the robustness of this kind of likelihood assessment. These can range from observational parameters, such as stellar age or O/Fe abundance ratios, to choices of theoretical models. For example, P(CHZ₂) can be calculated from models with more sophisticated HZ prescriptions. P(CHZ₂) itself can be used as a prior for calculating the likelihood of a higher level of detectability, perhaps combining it with the probability of geological activity given a measured age and stellar radionuclide abundances. For a subset of potentially habitable exoplanets, we have provided a quick and efficient method for observers to directly compare the estimated detectability potential for exoplanet host stars. This systematic statistical analysis will help to further narrow down the potentially habitable exoplanet candidates that could be targeted for in-depth follow-up characterization with upcoming space telescopes like JWST and future, more sensitive missions.