Beyond the Drake Equation: A Time-dependent Inventory of Habitable Planets and Life-bearing Worlds in the Solar Neighborhood

The search for habitable exoplanets and extraterrestrial life beyond the solar system is a topic of central interest for modern science and one of the most compelling and consequential endeavors for humankind. Simple life emerged on Earth within the first billion years of its habitable window, and the high frequency of terrestrial planets in the habitable zones (HZs) around GK dwarf stars inferred from NASA's Kepler observations (Bryson et al. 2021) invites the question of how often (if at all) life may have arisen on other worlds in the past. This pursuit will ultimately require statistical analyses of the population of habitable systems, in-depth studies of the climates of individual planets, and searches for chemical biomarkers (Schwieterman et al. 2018), and has motivated the development of the next generation of large ground-based facilities and instrumentation. The yield and characterization of Earth-like planets will be a primary science metric for future space-based flagship missions, but the optimal observational strategy for addressing the origin and properties of planetary systems and the prevalence of habitable exoplanets and life beyond the solar system remains unclear (e.g., Bean et al. 2017; Tasker et al. 2017; Sandora & Silk 2020; Truitt et al. 2020; Checlair et al. 2021; Sarkar 2022; Batalha et al. 2023). The gathering of comprehensive data for each individual system is impractical if not impossible, so a statistical perspective is necessary to prioritize targets for follow-up observations. In particular, one would like to identify—given a model of habitability and biosignature genesis—how potential biosignature yields change during the evolution of a stellar system as a function of stellar properties like age, mass, and metallicity.

This paper aims to present a theoretical framework for exoplanet population and astrobiology studies that may provide a better statistical understanding of the formation history, frequency, age, and metallicity distributions of different planet types around stars of different properties. Our approach may also establish a useful basis for testing hypotheses about habitable environments and life beyond the solar system, for gaining a sense of biosignature yields, and for informing future observational efforts and experiments. A well-known parameterization of the present-day abundance of life-bearing worlds in the Galaxy, N_ℓ , is represented by the first four terms in the probabilistic Drake equation (Drake 1965), which can be rewritten as

$$\begin{eqnarray}&&{N}_{{\ell }}={N}_{\mathrm{MS}}({t}_{0}){f}_{p}\,{n}_{e}\,{f}_{{\ell }}.\end{eqnarray} \tag{ 1 }$$

Here, N_MS(t₀) is the total number of stars that are on the main sequence today and can provide their planets a stable HZ, f_p is the fraction of these stars that have planetary systems, n_e is the average number per planetary system of Earth-size planets that are in the HZ, and f_ℓ is the subset of these rocky exoplanets that are "Earth-like" in a more detailed biochemical and geophysical sense and where simple life eventually arises. The first three terms (N_MS, f_p , and n_e ) in Equation (1) have already experimental measurements, and the fourth (f_ℓ ) is a conditional probability that may potentially be observable in the coming decades via spectroscopic searches for biosignature gases in exoplanet atmospheres (Schwieterman et al. 2018; Seager2018). In Drake's famous formulation, in order to estimate the number of active, communicative extraterrestrial civilizations around us today, the right-hand side of Equation (1) gets multiplied by the fraction f_i of life-bearing planets on which intelligent life emerges, times the percentage f_c of such civilizations that produces a detectable signal, times the fractional longevity f_L of a technological species.

The Drake equation and its "biosignature" version (Seager2018) amount to a pedagogical and organizational summary of the factors that may affect the likelihood of detecting technologically advanced civilizations or just simple microbial life evolving on a habitable planet. They are only meant to guide the observational inputs needed to make an educated guess, rather than provide a time-dependent mapping between star formation, environment, exoplanets, and life-harboring worlds. Their inherent limitations include both the lack of temporal structure (e.g., Ćirković 2004; Forgan 2009; Cai et al. 2021; Kipping 2021)—an assumption of uniformity with time that precludes the inclusion of evolutionary effects associated with, e.g., the star formation and chemical enrichment history of the local Galactic disk and the timeline of life emergence—as well as the difficulty of casting in a probabilistic argument the variety of phenomena and associated timescales that may influence anything quantified by probability f-factors and multiplicity n_e . Recent rapid developments in astrophysics and planetary sciences warrant a more informative and modern evolutionary framework, i.e., a rate-equation approach based on a system of first-order differential equations. These describe the changing rates of star, metal, planet, and habitable world formation over the history of a given stellar system, and can easily be adapted to incorporate the hierarchy of astrophysical and biological processes that regulate the age-dependent inventory of any key planet population.

The field of Galactic habitability and the formation history of Earth-like and giant planets in the Milky Way and the Universe as a whole have been research topics for more than two decades (e.g., Gonzalez et al. 2001; Lineweaver et al. 2004; Gowanlock et al. 2011; Behroozi & Peeples 2015; Gobat & Hong 2016; Zackrisson et al. 2016; Forgan et al. 2017; Balbi et al. 2020). Our approach expands upon some of the ideas employed in these early papers and develops new ones, focusing instead on the time-varying incidence of exoplanets and potentially habitable worlds over the chemo-population history of the solar neighborhood, the target of current and next-generation stellar and planetary surveys. This is the locale where more detailed calculations are justified by an avalanche of new data and actually needed in order to estimate, given a model of habitability and biosignature genesis, the relative biosignature yields among potential target stars.

The plan of this paper is as follows. In Section 2 we present the basic rationale and main ingredients of our modeling: the star formation history (SFH) and metallicity distribution function (MDF) in the solar vicinity, the planet occurrence rate around GK stars, and various metallicity-dependent effects. In Section 3 we cast and integrate our rate equations for the time evolution of the local abundance of dwarf stars and the giant planets and rocky planets in the HZ around them. We track giant planets to gauge the impact of their enhanced occurrence rate at higher host star metallicities relative to the weaker frequency–metallicity correlation of terrestrial planets. In Section 4 we extend our formalism to speculate about the formation history of life-harboring environments in the local volume under the hypothesis of a rapid abiogenesis process on Earth-like planets and estimate the prevalence of nearby biospheres in terms of the exoplanet census as a whole. Finally, we summarize our findings and conclusions in Section 5.

In order to provide absolute number counts in the solar vicinity we shall use the recent tally of main-sequence stars, giants, and white dwarfs within 100 pc of the Sun from the Gaia Early Data Release 3 (EDR3). The Gaia Catalog of Nearby Stars contains

$$\begin{eqnarray}&&{N}_{\star }({t}_{0})=\mathrm{331,312},\end{eqnarray} \tag{ 2 }$$

objects and is estimated to be >92% complete down to faint stellar type M9 (Gaia Collaboration et al. 2021). Apart from a minor correction (a fraction of a percent for the initial mass function (IMF) in Equation (4)) associated with the contribution of remnant neutron stars and black holes, N_⋆(t₀) represents the total number of stars ever formed in the solar neighborhood. Below, we shall use this normalization together with an SFH and an IMF to compute the number of main-sequence stars as a function of time.

2.1. SFH

Let ϕ(m) and ψ(t) be the (universal) IMF and SFH by number, respectively. The IMF and SFH are normalized so that ${\int }_{0}^{{t}_{0}}\psi (t){dt}=1={\int }_{{m}_{l}}^{{m}_{u}}\phi (m){dm}$. The number of stars that are on the main sequence at time t, N_MS(t), evolves at the rate

$$\begin{eqnarray}&&{\dot{N}}_{\mathrm{MS}}(t)={N}_{\star }({t}_{0})\int \phi (m)[\psi (t)-\psi (t-{t}_{\mathrm{MS}})]{dm},\end{eqnarray} \tag{ 3 }$$

where the dot denotes the time derivative, t_MS(m) < t is the main-sequence lifetime, and the second term in the square brackets corrects the rate of newly formed main-sequence stars for the number of stars that have evolved off the main sequence. In the following, we shall assume a Kroupa (2001) IMF

$$\begin{eqnarray}\begin{array}{c}\phi (m)=\left\{\begin{array}{cc}0.2530{m}^{-1.3} & \,\,\,(0.08\leqslant m\lt 0.5),\\ 0.1265{m}^{-2.3} & \,\,\,(0.5\leqslant m\lt 100),\end{array}\right.\end{array}\end{eqnarray} \tag{ 4 }$$

where m is measured in solar masses. The main-sequence lifetime t_MS can be computed using the analytical fitting formulae of Hurley et al. (2000) (based on the evolutionary tracks of Pols et al. 1998) as a function of m and metallicity Z (see Figure 3). ³ A G-type m = 1 main-sequence star, for example, has a lifetime of t_MS = 11 Gyr at solar metallicity, and t_MS = 6.5 Gyr at Z = 0.1 Z_⊙.

Zoom In Zoom Out Reset image size

The SFH of the solar neighborhood has been recently reconstructed by Alzate et al. (2021) using all 120,452 stars brighter than G = 15 mag within 100 pc of the Sun in the Gaia DR2 catalog. In broad agreement with previous determinations based on different techniques and data sets (e.g., Snaith et al. 2015; Mor et al. 2019; Ruiz-Lara et al. 2020), their results show two main early episodes of star formation: (1) a peak of activity occurring 10 Gyr ago that produced a significant number of stars with subsolar metallicities, followed by a star formation minimum (quenching) around 8 Gyr ago; and (2) a more recent burst about 5.5 Gyr ago. Since then, star formation has been declining until recent times, making stars with supersolar metallicities in short-lived bursts of activity.

Most low-metallicity 10 Gyr old stars belong to the thick Galactic disk population (e.g., Robin et al. 2014; Haywood et al. 2019), as opposed to the thin disk for the rest of them. Clear evidence of a thick-disk peak at age 9.8 ± 0.3 Gyr is also seen in the local white dwarf luminosity function by Fantin et al. (2019). An intense phase of star formation between 9 and 13 Gyr ago during the emergence of the thick disk, producing about as much mass in stars as that manufactured in the next 8 Gyr, and followed by a minimum in the star formation rate at an age of ∼8 Gyr, was also apparent in the SFH reconstruction of Snaith et al. (2014, 2015). In a 2 kpc bubble around the Sun, Ruiz-Lara et al. (2020) inferred three recent episodes of enhanced star formation dated ∼5.7, 1.9, and 1 Gyr, in synchrony with the estimated Sagittarius dwarf galaxy pericenter passages.

Figure 2 shows the marginalized posterior age distribution of solar neighborhood stars from Alzate et al. (2021) together with a reasonable reconstruction involving four Gaussians centered at ages of 10, 5.5, 2.0, and 0.7 Gyr, and having widths of 1.2, 0.9, 0.35, and 0.3 Gyr, respectively. In this reconstruction, which we use in the rest of this paper, about one-third of all stars belong to the old >9 Gyr population, and only 17% to the youngest, <3 Gyr component. There are periods of very little star formation around 7–8 Gyr ago and then again 3 Gyr ago. Note that, although the analysis by Alzate et al. (2021) uses a local sample, radial migration predicts that stars in the close solar vicinity may represent a mixture of stars born at various Galactocentric distances over the disk (see, e.g., Lian et al. 2022 and references therein).

Zoom In Zoom Out Reset image size

2.2. MDF and Age–Metallicity Relation

To study the influence of stellar metallicity on the occurrence rates of planets and planetary systems, we shall adopt here the MDF from the GALAH+TGAS spectroscopic survey of dwarf stars in the solar galactic zone (Buder et al. 2019). Figure 3 shows the data histogram and the corresponding best-fit skewed distribution, G(M), with moments 〈M〉 = −0.07, σ_M = 0.23, and Skew = −0.51, that accounts for the asymmetry of the extended metal-poor tail as well as the sharper truncation of the MDF on the metal-rich side. Here and below, we use the symbol M interchangeably with ${\mathrm{log}}_{10}(Z/{Z}_{\odot })$ for compactness, and the metallicity distribution g(Z) is related to the MDF in bins of M as $g(Z)=G(M)/({10}^{M}\mathrm{ln}10)$.

Zoom In Zoom Out Reset image size

Solar neighborhood stars exhibit an age–metallicity relation, such that young age correlates with high metallicity, a temporal sequence that is the fossil record of the enrichment history of the Galactic disk (e.g., Haywood et al. 2013; Hayden et al. 2015; Haywood et al. 2019). Observations have shown that this relation has significant scatter, attributed to the effects of radial migration and chemical mixing. Since our aim here is to characterize only the prevalent metallicity in each star formation episode correctly, we shall impose an age–metallicity relationship ignoring the dispersion around the mean—this is found to increase steadily with stellar age from 0.17 dex at age 2 Gyr to 0.35 dex at 13 Gyr (Buder et al. 2019). Within this framework, the fraction of stars that formed between time t and t + dt, ψ(t)dt, is then equal to the fraction of stars with metallicity between Z(t) and Z(t) + dZ, g(Z)dZ. The typical stellar metallicity therefore evolves with time as

$$\begin{eqnarray}&&\dot{Z}(t)=\psi (t)/g(Z).\end{eqnarray} \tag{ 5 }$$

The derived age–metallicity relation Z(t) is depicted in Figure 4. In broad agreement with previous results, it exhibits a rapid increase in metallicity at early epochs with an enrichment timescale of about 4 Gyr, followed by a slow evolution around solar values at ages between 4 and 10 Gyr and an upward trend toward supersolar metallicities at more recent times (e.g., Haywood et al. 2013; Snaith et al. 2015; Sharma et al. 2021).

Zoom In Zoom Out Reset image size

We note in passing that the age–metallicity distribution may actually consist of two distinct populations, an old and a younger sequence corresponding to the formation of, respectively, the thick and the thin disk (Nissen et al. 2020). An analysis of the implications of this scenario is postponed to future work.

2.3. Planet Occurrence Frequency Around FGK Stars

Exoplanet statistics in the inner regions of FGK dwarf stars have been investigated using the large and homogeneous sample from the Kepler mission (Thompson et al. 2018). Kepler planets commonly reside in multiplanet systems, and the integrated occurrence rate (the average number of planets per star) for exoplanets with radii in the range 1–20 R_⊕ and orbital periods up to P = 400 days is (Zhu & Dong 2021)

$$\begin{eqnarray}&&{\eta }_{P}\,=\,1.23\pm 0.06.\end{eqnarray} \tag{ 6 }$$

Earth-size exoplanets in Earth-like orbits are not well probed by Kepler, and estimates of their frequencies are more uncertain. A recent analysis by Bryson et al. (2020; see also Burke et al. 2015) yields

$$\begin{eqnarray}&&{\eta }_{1}\,=\,{0.015}_{-0.0007}^{+0.011},\end{eqnarray} \tag{ 7 }$$

for the occurrence rate around GK dwarf stars of terrestrial planets within 20% of Earth's orbital period and radius.

The planet radius–orbital period parameter space defining η₁ is a subset of the larger parameter space for η_⊕, the occurrence rate of Earth-size rocky planets in the HZ (hereafter "temperate terrestrial planets" or TTPs for short), roughly defined as the region around a Sun-like star in which a rocky planet with an Earth-like atmospheric composition can sustain liquid water on its surface (Kasting et al. 1993). The basic requirement for surface liquid water is predicated on a subset of the minimum conditions needed for a simple, microbial biosphere. Defining η_⊕ as the occurrence rate of TTPs with radii between 0.5 and 1.5 R_⊕ and orbiting stars with effective temperatures between 4800 and 6300 K, Bryson et al. (2021) recently derived

$$\begin{eqnarray}&&{0.37}_{-0.21}^{+0.48}\lt {\eta }_{\oplus }\lt {0.60}_{-0.36}^{+0.90},\end{eqnarray} \tag{ 8 }$$

where the errors reflect 68% confidence intervals and the lower and upper bounds correspond to different completeness corrections. This occurrence rate uses the conservative HZ estimates from Kopparapu et al. (2014).

Below, we shall adopt a fiducial present-day occurrence rate of η_⊕ = 0.24. This is the standard value used in forecasting TTP yields from direct-imaging future flagship missions like HabEx and LUVOIR, and is based on the NASA ExoPAG SAG13 meta-analysis of Kepler data (Kopparapu et al. 2018). Note that, in the language of Drake's equation (Equation (1)), η_⊕ ≡ f_p n_e .

For comparison, the frequency of giant gaseous planets with radii > 4 R_⊕ and orbital periods P < 400 days is estimated by Zhu & Dong (2021) to be

$$\begin{eqnarray}&&{\eta }_{\mathrm{GP}}\,=\,0.16\pm 0.015.\end{eqnarray} \tag{ 9 }$$

2.4. Dependence of Planet Frequency on Stellar Metallicity

In the context of core-accretion planet formation theory, metal-rich protoplanetary disks have enhanced surface densities of solids, leading to the more efficient formation of the rocky cores of gas giant planets (Pollack et al. 1996; Ida & Lin 2004). Is is well established observationally that metal-rich stars are more likely to host close-in giant planets (e.g., Fischer & Valenti 2005; Petigura et al. 2018; Zhu 2019), with an occurrence rate enhancement as a function of metallicity of the form

$$\begin{eqnarray}&&f{(Z)={ \mathcal A }(Z/{Z}_{\odot })}^{2}.\end{eqnarray} \tag{ 10 }$$

For a sample of stars with metallicity distribution g(Z), the normalization constant ${ \mathcal A }$ above is related to the integrated frequency of giant planets by Zhu et al. (2016)

$$\begin{eqnarray}&&{\eta }_{\mathrm{GP}}\,=\,\int g(Z)f(Z){dZ}.\end{eqnarray} \tag{ 11 }$$

With the adopted MDF (Figure 3), one derives ${ \mathcal A }\,=\,0.1369$.

Small planets show a weaker frequency–metallicity correlation (Sousa et al. 2008; Buchhave et al. 2012; Petigura et al. 2018). Most recently, Lu et al. (2020) were unable to confirm or reject a relationship between planet occurrence and host star metallicity for rocky planets with radii ≲ 2 R_⊕. In line with these analyses, no frequency–metallicity correlation for terrestrial planets will be assumed in our treatment. We shall also ignore possible correlations between small planet occurrence rates and α-element enhancement (Bashi et al. 2020).

2.5. Critical Metallicity for Terrestrial Planet Formation

In the core-accretion model of planet formation heavy elements are necessary to form the dust grains and planetesimals that build planetary cores. Johnson & Li (2012) estimated a minimum metallicity Z_c for planet formation by comparing the timescale for dust grain growth and settling to that for protoplanetary disk photoevaporation. They found that, for an Earth-size planet to form, a disk of surface density Σ(r) must have a metallicity

$$\begin{eqnarray}&&Z\gtrsim 0.1\,{Z}_{\odot }{\left[\displaystyle \frac{{\rm{\Sigma }}(1\,\mathrm{au})}{{10}^{3}\,{\rm{g}}\,{\mathrm{cm}}^{-2}}\right]}^{-0.48}.\end{eqnarray} \tag{ 12 }$$

Given the observed MDF (Figure 3), the existence of a fiducial metallicity floor Z_c = 0.1 Z_⊙ for the formation of terrestrial planets will only impact a small fraction of the census population in the solar neighborhood. Nevertheless, in our equations below we shall formally multiply the integrated occurrence rate of TTPs by the Heaviside function θ(Z–Z_c ).

2.6. Effect of Metallicity on the HZ

Low-metallicity stars have higher luminosities L_ZAMS and higher effective temperatures T_eff than metal-rich stars of the same mass (e.g., Tout et al. 1996). An m = 1 zero-age main-sequence (ZAMS) star of metallicity Z = 0.1 Z_⊙, for example, is 80% brighter, 13% hotter, and has a larger HZ than its solar-metallicity counterpart. Conversely, an F-type star with m = 1.5 and Z = 0.1 Z_⊙ has a main-sequence lifetime of only 1.8 Gyr (versus 2.7 Gyr at solar metallicity), and is unlikely to be a good candidate for harboring continuously habitable planets.

We assess the impact of stellar metallicity on the habitability of terrestrial exoplanets around GK stars (see also Danchi & Lopez 2013; Valle et al. 2014; Truitt et al. 2015) adopting the conservative HZ estimates of Kopparapu et al. (2014), where the inner and outer edges of the continuously HZ are defined by the "runaway greenhouse" (atmosphere becomes opaque to outgoing thermal radiation owing to excess amounts of H₂O) and "maximum greenhouse" (increased reflectivity of a thick CO₂ atmosphere wins out over greenhouse effect) limits. The 1D, radiative–convective, cloud-free climate models of Kopparapu et al. (2013) provide critical values for the effective flux S_eff—the normalized value of solar constant required to maintain a given surface temperature—as a function of the effective temperature of the host star

$$\begin{eqnarray}&&{S}_{\mathrm{eff}}\,=\,{S}_{\mathrm{eff},\odot }+{{aT}}_{\star }+{{bT}}_{\star }^{2}+{{cT}}_{\star }^{3}+{{dT}}_{\star }^{4},\end{eqnarray} \tag{ 13 }$$

where T_⋆ = T_eff − 5760 K and the coefficients (a, b, c, and d) for the runaway greenhouse and maximum greenhouse limits are listed in Kopparapu et al. (2014). Stellar effective temperatures and luminosities on the ZAMS were computed as a function of metallicity following Tout et al. (1996). The corresponding HZ distances can be calculated using the relation

$$\begin{eqnarray}&&{d}_{\mathrm{ZAMS}}\,=\,{\left(\displaystyle \frac{{L}_{\mathrm{ZAMS}}/{L}_{\odot }}{{S}_{\mathrm{eff}}}\right)}^{0.5}\,\mathrm{au},\end{eqnarray} \tag{ 14 }$$

where L_ZAMS/L_⊙ is the luminosity of the star in solar units. ⁴ Figure 5 depicts the predicted variations in HZ boundaries as a function of stellar metallicity for an m = 1 star. The effective stellar fluxes at the inner and outer edges of the HZ increase at low Zs because the calculated planetary albedos become higher as the star's radiation is shifted toward the blue, a dependence that is stronger at the outer HZ boundary because of the importance of Rayleigh scattering in dense CO₂ atmospheres (Kasting et al. 2014). Note, e.g., how a 1 M_⊕ planet at 1 au from an m = 1 ZAMS G-class dwarf of metallicity 0.1 Z_⊙ is just outside the inner edge of its host's HZ, as this is evaluated at 1.05–1.83 au.

Zoom In Zoom Out Reset image size

Figure 5 shows how, at low metallicities, the HZ moves further away from the host star and is about 25% wider. These differences in HZ boundaries may be relevant for present and upcoming planet-finding surveys around low-metallicity stars (Dedrick et al. 2020). Because old age correlates with low stellar metallicity, one might expect η_⊕ to be larger at earlier times if other factors, like the efficiency of planetary formation and planetary spacing, were equal. Given the relatively flat period distribution, ${df}/d\mathrm{ln}p\propto {p}^{0.2}$, observed for long-period Kepler planets (e.g., Bryson et al. 2020), however, the impact of such a shift in HZ boundaries on the predicted occurrence rates of TTPs around GK dwarfs in the young Galaxy should be rather minor.

2.7. HZ Lifetime

The boundaries of a radiative HZ are not temporally or spatially static but "migrate" outward over the course of a star's main-sequence phase. The secular increase in stellar luminosity results in a runaway greenhouse event, which, in the case of Earth, will cause the cessation of habitable conditions and the likely extinction of our biosphere about 1.5–2 Gyr in the future (e.g., Goldblatt & Watson 2012; De Sousa Mello & Santos Friaça 2020), well before the Sun becomes a red giant. The transitory nature of the residence of a TTP within an HZ has strong astrobiological implications and can be described by the time a planet, located at a given distance from a star, spends within the HZ (Danchi & Lopez 2013; Rushby et al. 2013; Waltham 2017). Here, we have used the formulae and best-fitting coefficients of Hurley et al. (2000) for the time-dependent luminosities and radii of stars on the main sequence to derive estimates for the change in the HZ boundaries over time.

Figure 6 (top panel) shows the evolution of the inner edge of the HZ, computed as before from the prescriptions of Kopparapu et al. (2014), as a function of time τ = t/t_MS for GK dwarf stars of different masses and metallicities. The point on each curve marks the characteristic "HZ lifetime," t_HZ(m, Z), when a hypothetical planet formed in the center of the HZ at the ZAMS stage enters into the hot zone of the host star, undergoes a runaway greenhouse event, and becomes uninhabitable (Rushby et al. 2013). ⁵ Lower-mass stars, while characterized by longer main-sequence lifetimes, have proportionally smaller t_HZ/t_MS ratios than higher-mass stars, the result of their lower rates of stellar luminosity evolution (Rushby et al. 2013). At fixed mass, main-sequence lifetimes are shorter and the total luminosity change over the main sequence is larger for lower-metallicity stars (Truitt et al. 2015). In the bottom panel of Figure 6 we compare the "HZ lifetime" with the main-sequence timescale as a function of stellar mass and metallicity. Over the plotted mass interval, the ratio t_HZ/t_MS ranges from 0.75 to 0.90 at solar metallicities, and from 0.65 to 0.75 at 0.1 Z_⊙. In absolute value terms, we find HZ lifetimes that are longer than the age of the Galaxy for m < 0.9 (Z = Z_⊙) and m < 0.75 (Z = 0.1 Z_⊙).

Zoom In Zoom Out Reset image size

Below, we shall assume that the evolution of main-sequence stars—rather than biogeochemical processes—is the only factor controlling the collapse of the HZ zone and the reduction of the biosphere lifespan. As the HZ expands outward due to the effects of stellar evolution, any planets that were initially beyond the boundaries of the HZ—so-called "cold start" icy planets—could potentially become habitable at later times as the HZ reaches them. We shall neglect this possibility below as the delayed habitability of such globally glaciated exoplanets remains dubious (e.g., Yang et al. 2017).

We can now cast our model for the time evolution of the exoplanet population in the solar neighborhood into a set of rate equations that can then be integrated as a function of time. Let us focus on K-class dwarfs with masses between m₁ = 0.45 and m₂ = 0.80 and main-sequence and HZ lifetimes that are typically longer than the age of the Galaxy at the metallicities of interest here. K stars may be better candidates in the search for biosignatures than G dwarfs, as they are more abundant, evolve less quickly on the main sequence, and provide their planets a stable HZ (e.g., Cuntz & Guinan 2016; Tuchow & Wright 2020). They also offer a longer photochemical lifetime of methane in the presence of oxygen compared to G dwarfs and, being dimmer, provide a better planet–star contrast ratio in direct-imaging observations (Arney 2019). ⁶ For the assumed Kroupa IMF, the fraction of stars that are classified as K-type is

$$\begin{eqnarray}&&{ \mathcal F }({m}_{1},{m}_{2})\equiv {\int }_{0.45}^{0.80}\phi (m){dm}=0.14.\end{eqnarray} \tag{ 15 }$$

Note that if the lower bound on habitability corresponded to the spectral type K5 instead (m = 0.65 M_⊙), the factor ${ \mathcal F }$ in Equation (15) would decrease by a factor of 3.5, while the inclusion of M dwarfs below 0.45 M_⊙ would boost the same integral by a factor of 6. Of course, for the occurrence rate of TTPs, what matters is the product ${\eta }_{\oplus }\times { \mathcal F }$, so one could group some of the uncertainties related to the lower habitability bound into the factor η_⊕.

3.1. Rate Equations

We can now track the evolution of the mean number of K-type stars in the solar neighborhood, N_K (t), and the abundance of giant planets and TTPs around them, N_GP(t) and N_⊕(t), by numerically integrating over time the corresponding rates

$$\begin{eqnarray}\begin{array}{rcl}{\dot{N}}_{K}(t) & = & {N}_{\star }({t}_{0})\psi (t){ \mathcal F },\\ {\dot{N}}_{\mathrm{GP}}(t) & = & f[Z(t)]{\dot{N}}_{K}(t),\\ {\dot{N}}_{\oplus }(t) & = & {\eta }_{\oplus }\theta [Z(t)-{Z}_{c}]{\dot{N}}_{K}(t).\end{array}\end{eqnarray} \tag{ 16 }$$

Here, in the equation for ${\dot{N}}_{\oplus }(t)$, we have assumed that exoplanets enter the HZ immediately after formation, and interpreted the parameter η_⊕ as the occurrence rate of TTPs around K-class stars on the ZAMS. ⁷ The equations above must be supplemented by Equation (5) for the evolution of the stellar metallicity Z(t).

Figure 7 shows the episodic formation history of exoplanets that results from the integration of Equation (16) for η_GP = 0.16 and η_⊕ = 0.24. In our "solar vicinity" sphere of 100 pc radius, there are currently about 11,000(η_⊕/0.24) TTPs around K dwarfs. They have a median age of 6.2 Gyr and 77% of them are older than the solar system. The minimum metallicity threshold for Earth-size planet formation of Johnson & Li (2012) does not significantly affect these numbers as the vast majority of star formation has taken place at Z > 0.1 Z_⊙. By contrast, the f(Z) modulation of the giant planet occurrence rate results in later typical formation times and shifts their median age to 3.9 Gyr, with terrestrial planets vastly outnumbering giant planets at early times.

Zoom In Zoom Out Reset image size

The early formation of TTPs in the solar vicinity occurred largely in two major episodes of enhanced star formation, starting with the emergence of the thick disk about 11 Gyr ago and followed by a second event that lasted 3.5 Gyr, peaked 5.5 Gyr ago, and involved more than 40% of the total stellar counts today. The five-planet system Kepler-444, orbiting a metal-poor 11.2 ± 1.0 Gyr old star, shows that thick-disk stars were indeed the hosts of some of the oldest terrestrial planets (Campante et al. 2015). The duration and size of the second major star formation surge suggest an external agent, perhaps a merger with a gas-rich satellite galaxy (Mor et al. 2019) or the first passage of the Sagittarius dwarf galaxy through the Milky Way's disk (Ruiz-Lara et al. 2020). Consistently with the Principle of Mediocrity, the solar system formed near the peak of this second episode. Over the last 4 Gyr the abundance of TTPs around K stars has increased by +24%, while that of giant planets has actually doubled.

In the coming decades, advanced space- and ground-based observatories will allow an unprecedented opportunity to probe the atmospheres and surfaces of TTPs in the search for signs of life or habitability. The discovery of extraterrestrial life would be a landmark moment in the history of science, with implications that would reverberate throughout all of society. Much of the early history of life on Earth has been dominated by methanogenic microorganisms, and methane in anoxic, Archean-like atmospheres is one of the most promising exoplanet spectroscopic biosignatures (e.g., Kaltenegger 2017; Schwieterman et al. 2018; Thompson et al. 2022). In spite of the recent swift developments in astrophysics and planetary sciences described in the previous sections, however, the probability of abiogenesis on Earth-like planets is currently unknown, as unknown are the characteristic timescales over which biochemical complexity actually evolves. The rapid emergence of life in the history of Earth has been used to argue for a high abiogenesis rate (e.g., Lineweaver & Davis 2002; Kipping 2020; Whitmire 2023). A Bayesian framework may naturally account for the anthropic bias that, if the timescale for intelligence evolution is long, life's early start may simply be a prerequisite to our existence, rather than evidence for simple primitive life being common on Earth-like worlds (Spiegel & Turner 2012). By means of this methodology, the recent analysis by Kipping (2020) shows that a fast abiogenesis scenario is at least three times more likely than a slow one.

In preparation for the next generation of space- and ground-based instruments, it seems interesting to conjecture on today's prevalence and time-varying incidence of microbial life-harboring worlds in the solar neighborhood under the hypothesis of a rapid abiogenesis process. We follow previous work (Spiegel & Turner 2012; Scharf & Cronin 2016; Chen & Kipping 2018; Kipping 2020) and describe abiogenesis as a stochastic Poisson process defined by a (uniform) rate parameter λ_ℓ —the mean number of abiogenesis events occurring per Earth-like planet in a fixed time span, which we set equal to 1 Gyr. The probability of achieving at least one successful abiogenesis event over a time interval $t-t^{\prime}$ since a given planet first became habitable at $t^{\prime}$ is then given by

$$\begin{eqnarray}&&{P}_{{\ell }}(t-t^{\prime} )=1-{e}^{-{\lambda }_{{\ell }}(t-t^{\prime} )}.\end{eqnarray} \tag{ 18 }$$

The time-dependent probability that life emerges on a TTP can then be expressed as the product ${P}_{{\ell }}(t-t^{\prime} )P({\ell }| \mathrm{HZ})$, i.e., we distinguish here between the population of planets that are "temperate" (i.e., Earth-size rocky planets in the continuously HZ) and the subset that are actually "habitable," i.e., "Earth-like" in a more detailed biochemical and geophysical sense, and where simple life will eventually arise. The number of these "Earth analogs" that formed (and became habitable) between time $t^{\prime}$ and $t^{\prime} +{dt}^{\prime}$ and where life emerged by time t is ${{dN}}_{{\ell }}(t)=P({\ell }| \mathrm{HZ}){P}_{{\ell }}(t-t^{\prime} ){\dot{N}}_{\oplus }(t^{\prime} ){dt}^{\prime}$. Assuming a probability P(ℓ∣HZ) that is independent of time, we can then write the mean number of life-hosting worlds present at time t as the convolution integral

$$\begin{eqnarray}&&{N}_{{\ell }}(t)=P({\ell }| \mathrm{HZ}){\int }_{0}^{t}{\dot{N}}_{\oplus }(t^{\prime} )\left[1-{e}^{-{\lambda }_{{\ell }}(t-t^{\prime} )}\right]{dt}^{\prime} .\end{eqnarray} \tag{ 19 }$$

Their formation rate must be given by the time derivative of Equation (19), yielding

$$\begin{eqnarray}&&{\dot{N}}_{{\ell }}(t)=P({\ell }| \mathrm{HZ}){\int }_{0}^{t}{\dot{N}}_{\oplus }(t^{\prime} ){\lambda }_{{\ell }}\,{e}^{-{\lambda }_{{\ell }}(t-t^{\prime} )}{dt}^{\prime} .\end{eqnarray} \tag{ 20 }$$

Note that, in the formalism of Drake's Equation (1), P(ℓ∣HZ) = f_ℓ in the limit of fast abiogenesis.

The oldest, generally accepted evidence for life on Earth comes from observations of microbial sulfate reduction in the 3.48 Gyr Dresser Formation (e.g., Lepot 2020). The Earth formed 4.54 ± 0.05 Gyr ago (Dalrymple 2001), and mineralogical evidence from detrital zircons indicates that liquid oceans must have been present on Earth's surface 4.404 ± 0.008 Gyr ago (Wilde et al. 2001). The maximum plausible value for the time interval over which at least one successful abiogenesis event occurred on Earth is therefore ≃0.9–1 Gyr. This is conservatively long compared to the maximum-likelihood timescale for life to first appear after conditions became habitable, ∼190 Myr, inferred by Kipping (2020), and much larger than the uncertainty as to when Earth became suitable for life, justifying our approximation of starting the "habitability clock" at formation.

We have integrated Equations (19) and (20) above assuming λ_ℓ = 1 Gyr⁻¹, and plotted in Figure 8 the resulting differential and cumulative number counts, ${\dot{N}}_{{\ell }}$ and N_ℓ . For illustration, we have assumed in the figure a conversion factor P(ℓ∣HZ) = 1 between TTPs and life-hosting Earths. Naturally, life would be abundant (again, we are concerned here with the appearance of the earliest life forms, not of intelligent life) in the solar vicinity if abiogenesis was fast and early Earth-like conditions existed and were relatively common on other worlds for 1 Gyr or more. The closest life-harboring exoplanet would be only 20 pc away if simple life arose as soon as it did on Earth in just 1% of TTPs around K stars. Conversely, Earth would be the only life-hosting planet in the solar neighborhood if abiogenesis was successful in about 1-in-10,000 TTPs. If microbial life is abundant it is also old, as it would have emerged more than 8 Gyr ago in about one-third of all life-bearing planets today. Note how the convolution integral in Equation (20) tends to smooth out the oscillations in ${\dot{N}}_{{\ell }}$ compared to the star and planet formation rates depicted in Figure 7, and that the assumed abiogenesis characteristic timescale of 1/λ_ℓ = 1 Gyr shifts the median age of the ∼10,000P(ℓ∣HZ) extrasolar biospheres predicted in the solar neighborhood to 5.7 Gyr.

Zoom In Zoom Out Reset image size

4.1. The Emergence of Oxygenic Atmospheres

A critical issue in the search for extraterrestrial life is whether Earth-like conditions lead to ecosystems that progressively oxygenate their planet atmospheres roughly following Earth's oxygenation history. The first oxygenation of Earth's atmosphere due to the emergence of photosynthesizing cyanobacteria happened about halfway through Earth's history (Luo et al. 2016), but O₂ rose irreversibly to near present atmospheric levels only about between 800 and 550 Myr ago, during the NOE (Och & Shields-Zhou 2012; Lyons et al. 2021), an event accompanied by major biological upgrades. While the precise timing of the NOE remains the subject of debate, our findings inevitably invite the question of whether and how often, given a habitable environment and following a successful abiogenesis event, the conditions for the beginnings of complex life may have arisen on exoplanets in the solar neighborhood. We can then consider a second stochastic process, labeled "O" for "oxygenation" and defined by a rate parameter λ_O , which can proceed only once abiogenesis ("ℓ") is successful. The inverse of λ_O is the characteristic timescale it takes for the earliest forms of life to evolve and produce an NOE-like juncture. Consider again an Earth-like planet that formed at time $t^{\prime}$. The joint probability density that abiogenesis was first successful at time t'' and was followed by an oxygenation event at time t is then given by

$$\begin{eqnarray}&&{p}_{O}(t^{\prime\prime} -t^{\prime} ,t-t^{\prime\prime} )={\lambda }_{{\ell }}{\lambda }_{O}\,{e}^{-{\lambda }_{{\ell }}(t^{\prime\prime} -t^{\prime} )-{\lambda }_{O}(t-t^{\prime\prime} )}.\end{eqnarray} \tag{ 21 }$$

The formation rate of simple life-hosting planets undergoing an NOE can then be written as

$$\begin{eqnarray}&&{\dot{N}}_{O}(t)=P({\ell }| \mathrm{HZ}){\int }_{0}^{t}{dt}^{\prime\prime} {\int }_{0}^{t^{\prime\prime} }{dt}^{\prime} {\dot{N}}_{\oplus }(t^{\prime} ){p}_{O}.\end{eqnarray} \tag{ 22 }$$

We have integrated this equation assuming ${\lambda }_{O}^{-1}\,=\,3.9\,$ Gyr and plotted the results in Figure 8. Because of the considerable delay between planet formation and NOE, the ∼7500P(ℓ∣HZ) worlds with strong oxygenic atmospheric biosignatures predicted to exist in the solar neighborhood have a formation rate that peaked 5 Gyr ago and a median age comparable to that of the solar system.

The search for habitable exoplanets and extraterrestrial life beyond Earth is one of the greatest scientific endeavors of all time. The high frequency of terrestrial planets in the HZs around dwarf stars implied by Kepler observations makes it timely to develop and explore new tools—beyond the probabilistic Drake equation—for statistical exoplanet population and astrobiology studies that may help direct future mission designs and observational efforts. In particular, one would like to understand—given a model of habitability and biosignature genesis—the formation history of simple life-harboring environments in the local volume and identify how potential atmospheric biosignature yields change as a function of stellar properties like age, mass, and metallicity.

The approach we have described in this work is based on a system of simple ordinary differential equations, rewritten below for the convenience of the reader

$$\begin{eqnarray}\begin{array}{rcl}{\dot{N}}_{K}(t) & = & {N}_{\star }({t}_{0})\psi (t){ \mathcal F },\\ \dot{Z}(t) & = & \psi (t)/g(Z),\\ {\dot{N}}_{\mathrm{GP}}(t) & = & f(Z){\dot{N}}_{K},\\ {\dot{N}}_{\oplus }(t) & = & {\eta }_{\oplus }\theta (Z-{Z}_{c}){\dot{N}}_{K},\\ {\dot{N}}_{{\ell }}(t) & = & P({\ell }| \mathrm{HZ}){\displaystyle \int }_{0}^{t}{\dot{N}}_{\oplus }(t^{\prime} ){\lambda }_{{\ell }}\,{e}^{-{\lambda }_{{\ell }}(t-t^{\prime} )}{dt}^{\prime} ,\\ {\dot{N}}_{O}(t) & = & P({\ell }| \mathrm{HZ}){\displaystyle \int }_{0}^{t}{dt}^{\prime\prime} {\displaystyle \int }_{0}^{t^{\prime\prime} }{dt}^{\prime} {\dot{N}}_{\oplus }(t^{\prime} ){p}_{O},\end{array}\end{eqnarray} \tag{ 23 }$$

which track the evolution of the mean number N_K of K-type stars in the solar neighborhood, their metallicity Z = Z(t), and the abundances N_GP, N_⊕, N_ℓ , and N_O of giant planets, TTPs, life-harboring worlds, and planets with oxygen-rich atmospheres, respectively These rate equations provide a time-dependent mapping between star formation, metal enrichment, and the occurrence of potentially habitable exoplanets over the chemo-population history of the solar neighborhood, and presents a useful basis for testing hypotheses about Earth-like environments and life beyond the solar system. The new framework can be easily adapted to incorporate the hierarchy of astrophysical and biological processes that regulate the age-dependent inventory of any key planet population.

We have numerically integrated the equations above adopting the recent tally of nearby stars (N_⋆) and white dwarfs from Gaia EDR3 (Gaia Collaboration et al. 2021), the episodic SFH (ψ) of the solar neighborhood as reconstructed by Alzate et al. (2021) and Ruiz-Lara et al. (2020), the MDF (g) from the GALAH+TGAS spectroscopic survey of dwarf stars in the solar galactic zone (Buder et al. 2019), and assuming an age–metallicity relation. In our model, the function f(Z) describes the metallicity modulation of the occurrence rate of giant gaseous planets (with integrated frequency today η_GP = 0.16; Zhu & Dong 2021), and we take η_⊕ = 0.24 for the fiducial occurrence rate of TTPs (Kopparapu et al. 2018) around K-class stars on the ZAMS. There is a minimum metallicity threshold for Earth-size planet formation of Z_c = 0.1 Z_⊙ (Johnson & Li 2012). Following earlier work (e.g., Spiegel & Turner 2012), we describe abiogenesis as a stochastic Poisson process defined by a (uniform) rate parameter λ_ℓ , and denote with P(ℓ∣HZ) the (constant) probability that a seemingly potentially habitable planet in the HZ was at early times "Earth-like" in a more detailed biochemical and geophysical sense and eventually became inhabited by life. A second stochastic process, an oxygenation event defined by a rate parameter λ_O , can proceed only once abiogenesis is successful.

Our main results can be summarized as follows.

• 1.  
  The formation of exoplanets in the solar vicinity followed two major events of enhanced star formation, starting with the emergence of the thick disk about 11 Gyr ago and followed by a second event that peaked 5.5 Gyr ago, lasted 3.5 Gyr, and produced more than 40% of all stars today. The solar system formed in the second star formation surge and was likely triggered by an external agent, perhaps a merger with a gas-rich satellite galaxy (Mor et al. 2019) or the first passage of the Sagittarius dwarf galaxy through the Milky Way's disk (Ruiz-Lara et al. 2020).
• 2.  
  Within 100 pc from the Sun, there are as many as 11,000(η_⊕/0.24) TTPs around K-type stars. About 77% of all TTPs in the solar neighborhood are older than the solar system.
• 3.  
  The metallicity modulation of the giant planet occurrence rate results in a later typical formation time, with TTPs vastly outnumbering giant planets at early times. Over the last 4 Gyr, the abundance of TTPs around K stars has increased by only +24%, while that of giant planets has doubled. The existence of a fiducial metallicity floor for the formation of terrestrial planets impacts only a small fraction of the census population in the solar neighborhood, as the vast majority of star formation has taken place at Z > 0.1 Z_⊙.
• 4.  
  The closest life-harboring Earth analog would be less than 20 pc away if microbial life arose as soon as it did on Earth in ≳1% of the TTPs around K stars. Conversely, Earth would be the only life-hosting planet in the solar neighborhood if abiogenesis was successful in about 1-in-10,000 TTPs. If simple life is abundant (fast abiogenesis with characteristic timescales 1/λ_ℓ = 1 Gyr), it is also old, as it would have emerged more than 8 Gyr ago in about one-third of all life-bearing planets today.

We finally note that errors in the number counts of exoplanets are likely dominated by planet occurrence rates (η_GP and η_⊕), which are uncertain at the ∼0.1–0.5 dex level due to incompleteness in long-period candidates. Comparable systematic errors may be associated with uncertainties in the IMF, SFH, and metallicity effects. Needless to say, in the case of life-hosting worlds, the precise values of P(ℓ∣HZ), λ_ℓ , and λ_O are currently unknown and remain a matter of speculation. Our work says nothing about how difficult or easy abiogenesis really is, a question that must ultimately be answered empirically. Still, given a model of habitability and biosignature genesis, our approach may provide a blueprint for assessing the prevalence of exoplanets and microbial life-harboring worlds over the chemo-population history of the solar neighborhood, gaining a sense of the atmospheric biosignature yields among potential target stars of different masses, ages, and metallicities, and guiding future observational efforts and experiments.

Support for this work was provided by NASA through grant 80NSSC21K027. We acknowledge useful discussions on this project with J. Alzate, M. Bisazza, F. Haardt, D. Lin, and R. Murray-Clay, and the hospitality of New York University Abu Dhabi during the completion of this study. The author would also like to thank the referee for a number of constructive comments and suggestions that greatly improved the paper.