Eta-Earth Revisited II: Deriving a Maximum Number of Earth-Like Habitats in the Galactic Disk

3.2.1. Uncertainties in the magnitude of the involved variables

Many factors that feed into the emergence and evolution of EHs cannot be scientifically quantified at present (see specifically Appendix 5), whereas others are only quantifiable up to a certain, low extent. Although our approach allows to set most of them equal to 1, several implemented parameters remain that are still highly uncertain. A prominent example is “Eta-Earth” itself, the occurrence rate of rocky exoplanets in the HZ of solar-like stars (see Sections 6.1.1 and B); an important variable through which we estimate the occurrence rate of rocky exoplanets in the HZCL. Literature values for η ⊕ since 2019 range over more than an order of magnitude (see Table 12), that is, from as low as η ⊕ ∼ 0.01 (Neil and Rogers, 2020) up to as high as η ⊕ ∼ 0.3 (Bryson et al., 2021), thereby still neglecting the specific error bars of these studies. This discrepancy partially stems from a relatively poor definition with no fixed lower and upper boundaries for orbital periods and planetary radii, indicating that scaling to the same boundary conditions slightly reduces this margin (see Kunimoto and Matthews, 2020, and Table 12).

However, further difficulties come into play here. The most crucial one relates to the fact that no rocky exoplanets were yet discovered in the HZ of solar-like stars, implying that observations from either lower mass stars or tighter orbital periods must be extrapolated. This is nontrivial since additional effects such as dependencies with stellar mass or temperature must be accounted for. Another often neglected effect that is increasingly taken into account in recent years relates to atmospheric erosion of primordial atmospheres on close-in orbits (Neil and Rogers, 2020; Pascucci et al., 2019), an important process that can lead to substantial overestimates of η ⊕ due to accidentally factoring in so-called evaporated cores. These are planets that were observed to be rocky close to the star but will be dominated by H₂-He-atmospheres in the HZ (Lopez and Rice, 2018; Pascucci et al., 2019).

Scaling from a specific range of orbital periods and radii to different parameter ranges, or from a specific stellar mass to another, as we need to perform for deriving the occurrence rate of rocky exoplanets in the HZCL, brings in further uncertainties, although these will be less pronounced than the inherent uncertainties in Eta-Earth itself. However, it will take quite some time until new instrumentation will discover rocky exoplanets in the HZs of solar-like stars and it will even take longer until robust occurrence rate statistics will exist on η ⊕ . Until then, we have to accept and implement relatively large uncertainties but need to highlight them properly.

Certainly, Eta-Earth is not the only parameter with high uncertainties. While the stellar parameters in our model are comparably well defined, any planetary parameters must be taken with caution. This is also true for the occurrence rate of planets with an appropriate water mass fraction and the frequency of large moons as both are mostly based on theoretical models but not on direct observations.8 Also, atmosphere and GHZ models are mostly based on theoretical considerations and on inherent assumptions that can lead to high uncertainties. As an example, are all relevant cooling agents implemented in an atmospheric model or are important molecules omitted, forgotten, or even unknown that may lead to a change in thermal stability (Section A)? What actually quantifies as a SN “sterilization” event is another relevant example.

These uncertainties are one of the reasons why we include literature reviews on various related topics (as mostly found in the appendices) and implement each of the parameters from comprehensive literature-based minimum to maximum values. They also highlight the possibility that our estimates can significantly change in the future when new models and observations are available.

3.2.2. Uncertainties in the choice of requirements themselves

Our derived formulation relies on the implementation of clusters of necessary requirements. However, whether a “necessary requirement” is indeed a necessary requirement is a category of caveat in itself (Ćirković, 2012). Although some of them are well accepted, others are not. The existence of a host star around which the rocky exoplanet orbits in the HZ is certainly well accepted, at least for the type of habitat we are considering, that is, for an “EH.” We define this as a rocky exoplanet within the habitable zone of complex life (HZCL) on which Earth-like N₂-O₂-dominated atmospheres with minor amounts of CO₂ can exist. Although this sounds almost tautological, this is not necessarily the case for other potential habitats such as the recently suggested Hycean worlds (Madhusudhan et al., 2021) for which no host star would be needed,9 for subsurface ocean worlds (Nimmo and Pappalardo, 2016) and planets with deep biospheres (Lingam and Loeb, 2020; McMahon et al., 2013), both of which do not have to orbit in the HZCL at all, or for brown dwarf (BD) habitats, which neither need a host star nor a conventional HZ, a rocky surface or an N₂-dominated atmosphere (Lingam and Loeb, 2019a). It is, however, less simple for other requirements.

Another example relates to the “right amount of water,” β pc H 2 O . It is well accepted that life as we know it (see discussion below) needs H₂O to evolve and thrive (Westall and Brack, 2018), even though this can already be quite subtle for certain extremophiles (Merino et al., 2019). However, it is certainly less accepted whether it is additionally required to have both subaerial land and oceans at the same time for “complex life” to evolve (another issue discussed below). In the appended Section 3.1, we give several plausible arguments for this assumption, but it is certainly not a well-established fact. Even if all the listed arguments in this section turn out to be true, there will still be the potential to resolve all of them by other means. For instance, if abiogenesis necessarily needs subaerial land (see appended Section 3.1), panspermia could in principle still transport life to such a planet (see Lingam et al., 2022a, and Section E.2.3). Our knowledge to provide definitive conclusions is yet certainly too poor. This in turn also relates to the first category described above: in case that complex life can commonly evolve on water worlds, our parameter range as based only on planets that have subaerial land and oceans may be too narrow.

There is another type of category, that is, the ones that are strongly debated and may finally not turn out to be a requirement at all. Most of them we did not implement into our model (some of these are discussed in Appendix 5), but one specific example we did, that is, the potential requirement of possessing a large moon. The initial argument in favor of the “Rare Moon” relates to the need for a large satellite to stabilize obliquity variations (Laskar et al., 1993). As it turned out, such an argument was not confirmed by simulations afterward (Waltham, 2019a), thereby illustrating that specific arguments must not necessarily be true. However, even though several further arguments in favor of a large moon exist (see appended Section 4.1), one could apply a similar reasoning as discussed above related to the water issue. Although we implemented a relatively “optimistic” parameter range for the occurrence rate of large moons (see Sections 6.3.1 and D.2), it may still overestimate its relevance if the requirement turns out to be negligible. On the contrary, if its importance will be confirmed in the future and if large moons are indeed rare, then our implemented parameter range could be too optimistic. Future observations can give important hints on this issue and on the uncertainties in the various parameters in general.

3.2.3. Methodological issues

Our formulation closely relates to the famous Drake equation (Drake, 1965) and therefore inherits some of its difficulties to a certain extent. One crucial problem with this type of equation is its lack of temporal structure (Ćirković, 2004; Lingam and Loeb, 2021). In its usual version, the Drake equation assumes a uniformity in time. As an example, the SFR of the galaxy feeds into the Drake equation via R ⋆ , the mean rate of SFR in the galaxy. However, the galactic SFR is not temporally uniform, as it was significantly higher in the past (see Snaith et al., 2015; Xiang et al., 2018, and Section 4.2.3). Assuming a mean rate could therefore lead to subsequent misconceptions. This is even more evident, if one takes the chemical “evolution” of the galaxy, as pointed out by Ćirković (2004). The metallicity in the galactic disk was much lower in the past making the formation of planets much more effective at present day (see Section 5.1.2). Averaging over the entire galactic history would therefore lead to wrong results. As a work around, one could divide the history of the galaxy into a stepwise function, one with too little metallicity in the past and another one with sufficient metallicity at present.

A similar assumption could be introduced for the GHZ by assuming that the Milky Way may have been uninhabitable in the past due to frequent SNs, again suggesting a stepwise function of habitability. In this case, f life , the fraction of planets on which life originates and evolves, could equal f life = 0 (or any other small number) for 0 < t ≤ t p , where t p signifies the phase transition from uninhabitable toward inhabitable with f life suddenly becoming a non-negligible number for the time interval t p < t ≤ t 0 , where t ₀ relates to the present day (Annis, 1999; Ćirković, 2004; Lingam and Loeb, 2021). These issues can lead to an under- or overestimate of habitable planets and life in today’s galaxy. It is therefore important to be accounted for properly.

However, our approach is not entirely equivalent to a uniform, step-less Drake equation since we effectively try to tackle this problem for specific parameters. As described in the following sections, we indeed implement the SFR, IMF, and stellar main-sequence lifetime into our model to arrive at a stellar age distribution for the Milky Way. This is of fundamental importance due to the aforementioned metallicity evolution and the thermal stability of N₂-dominated atmospheres. Because of the latter, we must account for stellar evolution, otherwise we could not calculate which stars allow for such atmospheres to be thermally stable (which is roughly favored for older stars; see appended Section 1). In the same manner, we also implement the evolution of galactic metallicity for considering the temporal distribution of rocky exoplanets. If we neglected these evolutionary parameters, we would overestimate the number of stars that can in principle host Earth-like atmospheres.

We do, however, not consider the evolution of SN rates, which indeed introduces a methodological uncertainty due to the spatial dynamics in the galactic disk. Drake equation-like formalisms cannot properly account for such dynamics, which is also true for our present framework. As an example, even though we only need to know the probability distribution of being sterilized by SNs for the present day, this does not assure us to choose the correct stellar systems for being sterilized because of neglecting stellar migration within the disk. A hypothetical star could, for example, migrate from a denser, more metal-rich region in the inner disk outward toward a metal-poorer region with far less detrimental SNs. This star could therefore have a sufficient metallicity for forming planets without being sterilized by SNs exactly because of its migration toward a quieter region. Since we ignore such migration, however, the same hypothetical star could have ended up being sterilized as it would have been if it had remained in the inner disk.

There are further evolutionary aspects we do not currently consider in our framework. One important example relates to the evolution of planetary atmospheres. Their composition (and pressure) can significantly change over geological timescales, thereby strongly affecting their thermal stability against atmospheric escape into space (see Sections 5.2.1 and A.2). Conceptually this is a crucial point in our formulation, although we do not consider it explicitly (but assume it implicitly). After the active phase of a star, when an N₂-dominated atmosphere theoretically becomes thermally stable at a planet, it must not necessarily be the case that such atmosphere will indeed evolve. We presently assume that any planet around a sufficiently weakly active star will host such an atmosphere—certainly a clear overestimate (see appended Section 5.2.1). We therefore set this requirement equal to unity, as we do with all requirements that are not implemented into our model.

Another crucial limitation of a Drake equation-like approach is the simple multiplication of parameters that may or may not be independent of each other (Ćirković, 2012). In reality, some of the parameters feeding into Relations 2, 3, and 4 can be positively or negatively correlated, implying that the outcome could be either an overestimate or underestimate, if known correlations are not addressed properly. We give some examples in Paper I, but here we emphasize an obvious positive correlation between two well accepted necessary requirements that we did indeed implement. These are the stellar metallicity and the occurrence rate, β HZCL , of rocky exoplanets in the HZCL. As we have knowledge about this positive correlation, we tackle it in our framework by weighting β HZCL for regions in the disk that show different metallicities than the solar neighborhood from which the occurrence rate is initially derived (see Section 5.1.2). Applying this methodology, as well as implementing evolutionary aspects such as the SFH and the evolution of metallicity, illustrates that our formalism allows for more complex calculations than simple multiplication of uniform values.

Potential correlations could in principle be tested and explored by performing statistical tests. This also includes the implementation of probability distribution functions instead of mostly taking point values or ranges (although these variables vary over space, time, and stellar mass in our model). This was already pointed out by several different authors who highlighted point variables instead of probability distribution functions to be another major limitation of the classical Drake equation (Glade et al., 2012; Lingam and Loeb, 2021; Maccone, 2010, 2011). However, implementing probability distribution functions is beyond the scope of the present study.

This brings us to a final methodological issue. As written above, we run three different cases for two types of atmospheres, a nominal case for which we always take mean values from the scientific literature, and minimum and maximum cases for which we always take minimum and maximum values, respectively. Therefore, both minimum and maximum cases are statistically highly unlikely since it is improbable that any chosen parameter will always reach either its theoretical minimum or maximum in reality. Based only on the implemented requirements, this would necessarily imply that the minimum cases are too low and the maximum cases are too high, with the actual values being closer toward the nominal cases. We further emphasize that this methodology still only results in a variation of the maximum number because several necessary requirements are not implemented into our framework but simply set equal to unity. Given this, it makes sense to not only consider maximum values for calculating the actual maximum number of EHs within the galaxy but to use the entire value range of each parameter from minimum to maximum for illustrating the uncertainties within the present scientific knowledge as well. If we were to only calculate the maximum number of EHs with strict maximum values from the scientific literature, the results would (i) pretend a strict theoretical maximum by neglecting any uncertainty range, and (ii) give unrealistically high maximum numbers as just described above.

3.2.4. Neglected factors

As already mentioned above, our approach ignores any potential requirements whose occurrence rates cannot be properly assessed with our current scientific knowledge. The occurrence rate of all these requirements is set equal to 1, which implies that our derived number of EHs must represent a maximum number by definition. The actual occurrence rate of EHs should therefore always be lower than the rate found by our formulation. An extensive list of potential galactic, stellar, planetary, and biological requirements ignored by our present study can be found in Appendix 5.

Our present formulation is further restricted in terms of potential habitats that it can presently cover, as it is designed for investigating the occurrence rate of EHs. It therefore does not cover any other hypothetical habitats such as subsurface ocean worlds, or more exotic ones such as the aforementioned Hycean worlds. However, this restriction is again by design. We are specifically interested in EHs because we can evaluate their existence at least partially based on some hard physical parameters. The most crucial of these parameters is the thermal stability of Earth-like atmospheres, that is, N₂-O₂-dominated atmospheres with a minor amount of CO₂. In this study, we emphasize the specific importance of the minor species CO₂, which not only relates to atmospheric stability but also to the toxicity limits of complex life (see below). Based on such relatively well understood atmospheric parameters, we can show that by far, not all stars are able to presently host EHs. This is a crucial point to emphasize because, based on specific known requirements, our formulation allows to derive some information on the prevalence of EHs. Even though this may seem similar to an anthropocentric ansatz, it is less clear how we can retrieve similar information on other types of habitats as our knowledge about them is much more restricted.

Another advantage of our approach relates to observability. Even though we are presently not at a stage where we can directly observe N₂-dominated atmospheres, we will be able to do so in the future. In the meantime, we can already start doing atmospheric statistics that will give important insights into the prevalence of EHs: If most planets in the HZCL have H₂- and/or CO₂-dominated atmospheres or even none at all, we can induce that N₂-dominated atmospheres will likely be rare. One can also infer that finding an N₂-O₂-dominated atmosphere will likely be an indication for life since these atmospheres function as a biosignature (see Paper I). This is an exciting prospect that makes deriving the maximum number of EHs in unison with their stellar/galactic distribution a worthwhile study to endeavor.

We should further explicitly highlight that we do not exclude the actual existence of other forms of habitats. If EHs are rare compared with the number of stars, it may well be that other forms of habitats could be more abundant. Although we are at present simply agnostic about them, we discuss a pathway for extending our formalism toward other forms of habitats in Section 7.8. We also point out that our formulation does not cover CETI like the Drake equation does. Also for this, we suggest a potential extension in Section 7.9.

Finally, we only investigate the galactic disk and neither the galactic bulge nor the galactic halo (as detailed in Section 4.1). We further neglect any space beyond our galaxy as it is not entirely trivial to expand our formalism toward a broader region of the Universe (Gonzalez, 2005). In addition, we neither calculate the evolution of habitability nor the maximum number of EHs throughout galactic history. Our calculations only relate to the present day. Finally, we want to emphasize that we do not consider parameters that may increase the maximum number of EHs. Obvious examples could be panspermia or the existence of hypothetically habitable exomoons. We do, however, discuss these issues in some detail in Section 7 and Appendix 5.

3.2.5. Some philosophical issues

As a last category, we briefly discuss certain philosophical issues related to our methodology, of which the definition of complex aerobic life is the most critical. We already discuss this issue in detail in Paper I and therefore keep it relatively short but it is of importance to reiterate it here. Our definition closely follows three studies that discuss certain atmospheric toxicity limits for complex life. These are Catling et al. (2005) related to O₂, Schwieterman et al. (2019a) related to CO₂, and Ramirez (2020) related to CO₂ and N₂, that is, the three prime atmospheric species for our framework.

Based on these studies, we use the terms “complex (aerobic) life,” “life as we know it,” “advanced metazoans,” and “animal life” synonymously to mean millimeter- to meter-sized carbon-based heterotrophs that are mobile and contain a blood-circulatory-like system comparable with advanced metazoans.10 Whenever we mean other forms of life in our study, we spell it out explicitly. We also note that focusing on habitats suitable to animal-like complex life, although it may seem arbitrary at first glance, may actually be well justified. On Earth, animals are the only complex multicellular eukaryotes capable of phagocytosis and animals are effective “ecosystem engineers.” For these reasons, Lingam and Loeb (2021) suggest that the evolution of animal multicellularity is one of the five, potentially universal, key critical steps for the emergence of technological intelligence on Earth and on exoplanets in general.11

Catling et al. (2005) precisely outline the universal importance of atmospheric O₂ for such complex life by showing that a certain O₂ partial pressure of pO 2 ≳ 0.1 PAL (with PAL as Earth’s present atmospheric level of O₂) is needed to allow for the emergence of complex life above a size of about 1 cm in diameter based on a crossover of the maximum physical size for O₂ diffusion and the minimum physical size for the emergence of a blood circulatory system. Based on such physical and energetic arguments, one can therefore expect that complex life needs similar O₂ partial pressures also on other worlds (Catling et al., 2005). Note, however, that the lower O₂ partial pressure threshold for the emergence of technological species may be even more restrictive than the one for complex aerobic life, since this likely needs pO 2 ≳ 0.86 PAL (i.e., a mixing ratio of x O 2 ≥ 18 % for a 1 bar atmosphere) based on certain combustion limits (Balbi and Frank, 2024). The CO₂ toxicity levels of advanced metazoans on Earth are further outlined in Schwieterman et al. (2019a) and Ramirez (2020) who both find the maximum level of pCO₂ to be about 0.1 bar (i.e., 10% CO₂ in a 1 bar atmosphere). Ramirez (2020) additionally investigates N₂ toxicity limits for complex life and finds maximum pressures of about pN 2 ∼ 2 bar to be lethal.

In this study, we need to highlight, however, that extraterrestrial life must not necessarily evolve toward similar toxicity limits but may find additional pathways to cope with toxic atmospheres, especially if the evolutionary pressure of survival is high. We nevertheless find these pressure limits to be a useful starting point for evaluating the prevalence of EHs in the galaxy even if such limits may vary substantially on other worlds. Until no other biospheres are found and investigated, it is reasonable to evaluate atmospheres within aforementioned limits as long as it is made clear that these can be substantially different for extraterrestrial complex life.

Whenever we talk about Earth-like atmospheres, we implicitly mean N₂-dominated atmospheres with O₂ as a second, less abundant main species, and with a minor amount of CO₂. Such atmospheres can certainly include other trace species and the exact ratios of the three main species are not fixed as long as the CO₂ mixing ratio is either below x CO 2 = 10% or x CO 2 = 1%. We have chosen these specific values because (i) x CO 2 = 10% CO₂ coincides with pCO 2 = 0.1 bar for 1 bar total pressure, that is, the aforementioned toxicity limit for complex life, (ii) HZCL boundaries exist for exactly these mixing ratios related to a 1 bar atmosphere, and (iii) atmospheric models of N₂-dominated atmospheres exist for the same mixing ratios as well. In addition, the x CO 2 = 1% CO₂ case more closely resembles the actual atmosphere of the Earth, which is interesting in itself. We further emphasize that 10% and 1% CO₂ mixing ratios can always be regarded as maximum values and we will henceforth use 10% and 1% CO₂ synonymously to 0.1 bar and 0.01 bar CO₂.

This specific definition of an atmosphere also relates to our definition of what an EH actually constitutes. As we have already stated, we define an EH as a rocky exoplanet within the HZCL on which Earth-like N₂-O₂-dominated atmospheres with minor amounts of CO₂ can exist. This definition in principle allows for the possible existence of complex aerobic life that obeys to similar limits as advanced metazoans here on Earth, but it is important to note that such putative organisms do not have to actually live on such a planet. However, the proposition that the Earth-like atmospheres themselves, that is, the simultaneous existence of N₂ and O₂ with minor amounts of CO₂, act as a biosignature makes EHs extremely relevant for astrobiology in itself, regardless of whether complex life actually exists there or not.

Although EHs are close to the definition of Class I habitats in Lammer et al. (2009),12 both are not entirely the same since the latter does not directly relate to atmospheric composition. It does, however, relate to certain stellar and geophysical conditions needed for such habitats to evolve, which can be regarded to be semantically equivalent with our meaning of necessary requirement. Class I habitats further exclude ocean worlds explicitly on which the ocean is not in direct contact with the mantle silicates due to the formation of a high-pressure (HP) ice layer between both (such planets are defined as Class IV habitats in Lammer et al., 2009). We also exclude them from our definition, but this does not relate to ocean planets without subaerial land on which an HP ice layer may not exist. We discuss arguments on whether these planets could constitute an EH in detail in appended Section 3.1. However, we also note that Lammer (2013) refined the definition of Class I habitats, so that (i) it does exclude ocean worlds without subaerial land entirely (with ocean worlds constituting Class V habitats) and (ii) N₂ as main species should be present in their atmospheres.13

If one takes the definition of EHs (or Class I habitats), it seems logical that it becomes narrower the more is scientifically known about said certain stellar and geophysical conditions. If we, for instance, could exclude planets without subaerial land based on future scientific knowledge, this will narrow down its definition and hence the number of planets that can evolve into an EH. One must, however, take care to not overspecify its definition, otherwise the definition of EH will be indistinguishable from Earth, and Earth will by definition be the only EH in existence. It is therefore important to keep a balance between a too narrow, tautological definition and a broader useful one that allows searching for planets on which complex aerobic life may evolve based on identifying and testing relevant and scientifically quantifiable requirements.

Another issue is related to the GHZ. It is important to highlight that the models we utilize for the sterilization rates of SNs were performed for complex life as defined above. As we discuss in Section 5.1, the effects of SNs would be less pronounced if we would only consider extremophiles or microbial life. We also want to mention that the definition of rocky exoplanet, and hence of “Earth-like” itself, seems to be slightly contentious, as there is no fixed mass or radius range that applies to the term rocky exoplanet. We discuss this in detail in Sections B and 6.1. Finally, we also emphasize that our results have certain philosophical implications such as for the Copernican Principle, the Fermi Paradox, and the Search for and Messaging to Extraterrestrial Life, that is, for SETI and METI. We discuss these potential implications in Section 7.

4.1.1. The various components of the Milky Way

The Milky Way (see reviews by Bland-Hawthorn and Gerhard, 2016; Helmi, 2020) is a typical disk galaxy and can be divided into different regions, that is:

4.1.2. Excluding the galactic halo

As written above, the galactic halo mainly consists of metal-poor and old stars. Zuo et al. (2017), for example, can reproduce the MDF in the galactic halo well by separating it into three distinct groups with peak metallicities of (Fe/H) approximately −0.63, −1.45, and −2.0, where the two metal-poorer components correspond to the inner and outer halo, respectively.14 The additional, relatively metal-rich component, which corresponds to substructures within the galactic halo such as clouds, and streams such as the Sagittarius stream (Belokurov et al., 2007; Grillmair and Carlin, 2016; Koposov et al., 2012), covers about 10% of halo stars (Zuo et al., 2017). Similar results were found by other surveys, for example, (Fe/H) approximately −1.2 ± 0.3, and −2.0 ± 0.2 by Liu et al. (2018) for the inner and outer galactic halos, respectively.

It seems unlikely that many of them will provide habitable conditions at present, however. As already pointed out by Lineweaver et al. (2004), many of the halo stars will have metallicities that are too low to host rocky exoplanets, a threshold that might be somewhere around [ Fe/H ] ∼ − 0.5 (Johnson and Li, 2012, see Section 5.1.2), a value that is roughly similar to the aforementioned metal-rich component of the galactic halo. If we assume a stellar mass of ≈ 1.3 × 10 9 M e (Deason et al., 2019; Mackereth et al., 2019) for the entire halo, the 10% component of metal-rich stars will equate to ≈ 1.3 × 10 8 M e . That number is only about ∼0.3% of the entire stellar disk mass, if we take the recently calculated value of by Xiang et al. (2018), which is around the average stellar mass of all disk models (see, Wang et al., 2020). Roughly 3 out of 1000 stars15 that could in principle host rocky exoplanets in the combined region of disk and halo will therefore belong to the halo. However, this is a relatively optimistic estimate, as the peak metallicity of (Fe/H), approximately −0.63 for the metal-rich halo component, still suggests that a high number of stars lie below (Fe/H) = −0.5.

In addition, the high majority of stars within the galactic halo, including the inner halo, are of an age of ∼10–12 Gyr (Jofré and Weiss, 2011) while the metal-rich stellar component of the Sagittarius stream holds ages of ∼9.5–11 Gyr (Carollo et al., 2016). Similar ages can be found for other metal-rich halo regions such as the Styx and Orphan streams, both of which were found to cover ages of ∼10–11 Gyr (Carollo et al., 2016). These stellar ages are significant, and one can expect that many of the rare planets that formed around these stars do not show geological activity at the present day, specifically if one considers that planets with an Earth-like radiogenic heat budget will be geologically active for about 6 Gyr (Mojzsis, 2021). These worlds would therefore provide relatively limited conditions for the evolution and maintenance of complex life. In fact, these stellar ages are even too old to generally allow inhabited planets around G-type stars based on their average main-sequence lifetime and bolometric luminosity evolution. EHs in the galactic halo must therefore logically be restricted either to old low-mass stars, which pose further difficulties (as discussed in later sections), or to one of the few younger ones that were either ejected out of the galactic disk (Faltová et al., 2023) or formed at younger ages under as of yet poorly constrained conditions (Bellazzini et al., 2019).

No study has explicitly modeled the habitability of the galactic halo, so no potential distribution can be incorporated within our model (as we do for the disk in the following subsections). However, even if there are stars within the halo that currently allow for the existence of an EH, their numbers will be minuscule, most likely significantly below 1 permille of the entire population, at least based on the considerations above. Due to their average distance, the potential to observe and characterize some of them in the near future may further be relatively elusive. Based on this reasoning, we therefore do not consider the galactic halo within this study.

4.1.3. Excluding the galactic bulge

The story looks a bit different for the galactic bulge. Depending on the specific study (Licquia and Newman, 2015; McMillan, 2017, 2011), its stellar mass amounts to about ∼20% of the entire galactic stellar mass. So, if we assume the bulge to be habitable, it cannot be neglected as long as we want to estimate the maximum number of EHs within the entire Milky Way. Whether the bulge is indeed habitable, however, is debated (Balbi et al., 2020), as outlined in the following.

The habitability of different galactic regions, which is often subsumed under the concept of the so-called GHZ, depends on various factors (Kaib, 2018) that are usually broadly divided into three (nonexhaustive) areas, that is,

high-energetic events such as AGN, SNs, and GRBs that both sterilize a planet and erode its atmosphere;

Several specific problems may arise for the bulge, particularly within (i) and (ii).

Since the stellar density in the bulge is significantly higher than in the disk and halo regions, it can be expected that a planet suffers from detrimental SNs more often than a planet in the other disk regions. As estimated by Gehrels et al. (2003), an SN within r ≤ 8 kpc can significantly affect a planet’s habitability by depleting about 30% of its ozone layer and thereby doubling the incident UV irradiation; a value that can trigger mass extinctions of land-based life (Melott and Thomas, 2011). While the occurrence rate of such an event at Earth is calculated to be ∼ 1.5 Gyr^–1 based on the average galactic rate of core-collapse SNs (Gehrels et al., 2003), Balbi et al. (2020) estimated that such extinction events may have a frequency of ∼ 40 − 110 Gyr^–1 at an average planet within the bulge. However, based on the resilience of Ecdysozoa, a group of protostome animals that include Milnesium tardigradum (i.e., the highly resilient tardigrades, also known by the name of water bears), Sloan et al. (2017) calculated that it needs an SN within ≤ 0.04 pc to completely sterilize a habitable planet. By considering such boundary condition, Balbi et al. (2020) estimated that the expected number of such events closer than 2 kpc from the center is significantly smaller, that is, ∼ 0.6 × 10 − 6 Gyr^–1 to ∼ 1.4 × 10 − 5 Gyr^–1. However, it should be further noted that tardigrades are exceptionally resilient survivors and are clearly nonrepresentative of the average animal. While the radiation threshold LD₅₀, at which 50% of the entire investigated population dies, is in the order of 5000 Gray for M. tardigradum within 48 h (Horikawa et al., 2006), this value is around a few Gray for typical animals such as rats, frogs, and deer (ICRP Publication 108, 2008), that is, ∼1000 times lower.

Apart from SNs, GRBs (Gowanlock, 2016; Melott et al., 2004; Piran and Jimenez, 2014; Scalo and Wheeler, 2002; Thomas et al., 2021; Thomas et al., 2005b, 2005a) could be another threat within the bulge. Piran and Jimenez (2014), for example, found that these could render the galactic center inhospitable to life. However, a study by Gowanlock (2016) that considers the SFH and metallicity evolution of the galaxy suggests a more complex picture. If GRBs correlate with low-metallicity environments (Jimenez and Piran, 2013), predominantly the metal-poor outskirts of the Milky Way will presently be affected by GRBs. If GRBs, on the contrary, exclusively follow the galactic SFH, most of them will occur in the galactic center.

The supermassive black hole (SMBH) Sagittarius A ⋆ (Sgr A ⋆ ) presents another potential threat to habitability (Ambrifi et al., 2022; Balbi et al., 2020; Balbi and Tombesi, 2017; Forbes and Loeb, 2018; Heinz, 2022; Lingam et al., 2019; Wisłocka et al., 2019). Even though its early phase as an AGN had a likely initial duration of ≲ 10 8 Myr (Marconi et al., 2004), other black hole (BH)-related processes may still present a threat to the habitability of the bulge. The frequency of so-called tidal disruption events, that is, the disruption of a star that passes too close to the SMBH, is estimated to be ≲ 100 Myr^–1 (Komossa, 2015; Stone et al., 2020). Pacetti et al. (2020) found that the induced radiation from such an event could significantly harm habitability within distances of ∼ 0.1 − 1.0 kpc from the galactic center. In addition, stochastic processes such as gravitational instabilities can frequently activate SMBHs in the mass range of Sgr A ⋆ every 10 7 − 10 8 years with an outburst duration of ∼ 10 5 years (Hopkins and Hernquist, 2006). Such outbursts could again render planets uninhabitable within ≲ 1 kpc (Amaro-Seoane and Chen, 2019; Ambrifi et al., 2022). There is indeed evidence that an outburst of X-ray radiation took place around 2 to 8 Myr ago (Bland-Hawthorn et al., 2013; Nayakshin and Cuadra, 2005), illustrating that Sgr A ⋆ cannot be neglected when studying the habitability of the bulge, not even at present day.

As a side note, it should be mentioned that the effect of the XUV flux from an SMBH on the atmospheric mass loss of a rocky exoplanet may be underestimated. Balbi and Tombesi (2017) and Balbi et al. (2020) calculate this process with the well-known energy-limited escape esquation (see Equation 25 in Section 6.2.1), in which the square of the photospheric radius, R ph , and the extreme UV (EUV) flux absorption radius, R XUV , enter the numerator16 (Watson et al., 1981). However, by simply assuming R ph 3 instead of R ph 2 R XUV , as many studies do, they do not consider the XUV flux-induced extreme expansion of a hydrodynamically expanding atmosphere. This leads to a substantial underestimate of this loss process, implying that the deleterious effects from an SMBH could be substantially higher than calculated by these studies (see also Appendix 1).

Because of the high stellar density in the galactic center, stellar systems are susceptible to orbital perturbations by nearby stars (Balbi et al., 2020; Bojnordi Arbab and Rahvar, 2021; de Juan Ovelar et al., 2012; Jiménez-Torres et al., 2013; McTier et al., 2020). McTier et al. (2020), for example, found that 8 out of 10 bulge stars experience stellar encounters within 1000 AU at a rate of ∼ 1 Gyr^–1, while about half of all bulge stars even experience >35 Gyr^–1. Such close encounters can either significantly perturb protoplanetary disks around such stars (Bhandare and Pfalzner, 2019; Cai et al., 2019; Li et al., 2019), or destabilize and even eject planets from their hosts (Vincke and Pfalzner, 2018; Winter et al., 2018b, 2018a). See appended Section 5.1.1 for a more detailed discussion.

Finally, we point out a potentially positive aspect of habitability related to the bulge. Apart from its generally high density of stars, which can be regarded as a positive aspect in itself at least related to the number of systems that can be examined per volume of space, the advantage for panspermia, that is, for transferring life from one stellar system to another, should particularly be highlighted (Adams and Spergel, 2005; Balbi et al., 2020; Belbruno et al., 2012; Chen et al., 2018; Gobat et al., 2021; Chen, 2021).17 In the bulge, the travel time of an ejecta for reaching the nearest star can be smaller by at least an order of magnitude compared with the galactic disk; it may hence be within the estimated survival time of certain extremophiles (Balbi et al., 2020; Ginsburg et al., 2018).

The arguments of the last paragraphs do not paint an entirely coherent picture or generally prove the complete noninhabitability of the galactic bulge, at least not for microbial life and/or extremophiles. However, they indicate that other galactic regions, particularly within the disk, might provide far more favorable habitable environments than the galactic center. Because of this and the lack of studies on the SFH and the distribution of SNs and metallicity within the bulge, we presently restrict ourselves to estimating the maximum number of EHs in the galactic disk.

To derive such a maximum number, however, we need to calculate the spatial distribution, age, and total number of presently existing stars in the galactic disk. To achieve such a distribution, we first need to implement the IMF and SFH into our model.

4.2.1. The IMF

The IMF describes the initial distribution of stellar masses, M ⋆ , for a population of stars, N ⋆ . It was first prominently introduced by Salpeter (1955) and can in its simplest form be written as a power law such thatd N ⋆ d M ⋆ = ξ 0 M ⋆ α ,(5) where ξ ₀ is a normalization constant and α is the power law index, which is α = − 2.35 in the case of the Salpeter-IMF (Salpeter, 1955). By integrating Equation 5, that is,∫ d N ⋆ = ξ 0 ∫ M ⋆ low M ⋆ up M ⋆ α d M ⋆ = ξ 0 α + 1 [ M ⋆ α + 1 ] M ⋆ low M ⋆ up ,(6) one will get the fraction of stars between the lower and upper stellar masses, M ⋆ low and M ⋆ up , respectively. For obtaining the distribution of the number of stars, δ N ⋆ , within certain mass bins, δ M ⋆ , we have to set N ⋆ = 1 and calculate ξ ₀ by inserting the minimum and maximum stellar masses for M ⋆ low and M ⋆ up .

The lower stellar mass limit, M ⋆ low , is usually constrained by the hydrogen burning limit to be around 0.07 and 0.08 M ⊙ (Kroupa, 2001; Kroupa et al., 2013) or set to be 0.1 M ⊙ (Chabrier, 2005, 2003; Salpeter, 1955).18 The upper stellar mass limit, M ⋆ up , usually resides between 100 M ⊙ and 150 M ⊙ (Kroupa et al., 2013), but its exact chosen value has no notable influence on the distribution because stars with higher mass are becoming increasingly rare. The lower limit, on the contrary, shifts the distribution toward smaller stellar masses, because low-mass stars are dominating the distribution (see also Table 1). Such a shift leads to a lower average stellar mass and therefore to more stars within the disk. However, since we only consider stellar masses from 0.1 M ⊙ to 1.25 M ⊙ (see Section 5.2), varying M ⋆ low between 0.07 M ⊙ and 0.1 M ⊙ does not substantially alter the final number of stars that we consider for deriving a maximum number of EHs (see also Table 6).19

The IMF was subsequently further investigated by several different authors (Chabrier, 2005, 2003; Kroupa et al., 1993, 2013; Kroupa, 2001; Miller and Scalo, 1979) to better fit the refined observational data. Kroupa et al. (2013) provide a thorough review of the IMF and propose a two-part power-law “canonical IMF” for masses between M ⋆ low = 0.07 M ⊙ and M ⋆ up = 150 M ⊙ based on the IMF by Kroupa (2001). We implement this canonical IMF in our nominal case because it is commonly used within various relevant models of the galactic disk mass or SFH. It can be written as (Kroupa et al., 2013)d N ⋆ d M ⋆ = ξ 0 { ( M ⋆ 0.07 ) − 1.3 ± 0.3 , 0.07 ≤ M ⋆ ≤ 0.5 M ⊙ [ ( 0.5 0.07 ) − 1.3 ± 0.3 ] ( M ⋆ 0.5 ) − 2.3 ± 0.36 , 0.5 ≤ M ⋆ ≤ 150 M ⊙ .(7)

By integrating this relation, we obtain the number of stars, δ N ⋆ , within each stellar mass bin, δ M ⋆ . This can be displayed as a cumulative mass fraction of stars, as plotted in Figure 1 for stars with stellar masses up to M ⋆ = 2.0 M ⊙ . There, the solid lines depict mass distributions based on various IMFs with M ⋆ low = 0.07 M ⊙ and M ⋆ up = 100 M ⊙ as boundary conditions, while the dashed and dotted lines show the same IMFs but for M ⋆ low = 0.08 M ⋆ e and M ⋆ low = 0.1 M ⊙ , respectively. As can be seen, the change in M ⋆ low mostly affects the distribution within the spectral class of M dwarfs (see also Table 1). For our nominal case we use M ⋆ low = 0.08 M ⊙ and M ⋆ up = 100 M ⊙ as suggested by Kroupa (2001), while we implement M ⋆ low = 0.07 M ⊙ and M ⋆ low = 0.1 M ⊙ for our minimum and maximum cases, respectively (with M ⋆ up always set to be 100 M ⊙ ). We further restrict ourselves to the IMFs by Kroupa (2001) for the nominal and minimum cases, and Chabrier (2003) for the maximum case. We do not consider the early, seminal IMF by Salpeter (1955) since it was already shown by Paresce and De Marchi (2000) that a single power-law mass function cannot reproduce the whole stellar distribution but overestimates stars below M ⋆ ∼ 0.5 M ⊙ .

OPEN IN VIEWER

FIG. 1.

Different cumulative stellar mass distributions for different initial mass functions (IMFs), that is, of an initial population of stars, as calculated through the empirical power laws by Salpeter (1955), Miller and Scalo (1979), Kroupa (2001)/ Kroupa et al. (2013), and Chabrier (2003), displayed for stellar masses up to 2.0 M ⊙ . The solid lines were calculated between stellar mass limits of 0.07 − 100 M ⊙ , while the dashed and dotted lines show calculations for 0.08 − 100 M ⊙ , and 0.1 − 100 M ⊙ , respectively. The dotted vertical lines give the borders of the different stellar spectral classes (after Habets and Heintze, 1981).

We should finally mention that the IMF may not be constant but varies over time. Li et al. (2023), for instance, found that the relative fraction of low-mass M dwarfs that form during an episode of star formation increases with stellar metallicity and hence galactic age. Based on correction factors by Liu (2019) and Moe et al. (2019), Li et al. (2023) also provide an IMF correction for binary systems that may further increase the relative fraction of young low-mass stars. This indicates that a higher number of M dwarfs form at present than at earlier ages of the galaxy. Since we implement the canonical IMF and do not consider such variability, this may slightly affect our results. If we implemented a variable IMF instead, our stellar distribution would contain a higher number of young M dwarfs, while the fraction of stars in the other spectral classes would be slightly lower. Such alteration would potentially lead to a slightly lower maximum number of EHs since a larger number of stars would fail to meet A at i , that is, the necessary requirement of permitting the long-term stability of an N₂-O₂-dominated atmosphere (see specifically Section 5.2).

Since the IMF describes the initial stellar distribution, we need to implement the stellar main-sequence lifetime as a function of stellar mass to calculate the potential survival of a star until the present day. For the main-sequence lifetime, we follow the same approach as Westby and Conselice (2020) who combine the estimated main-sequence lifetime of the Sun of about 10 Gyr (Schröder and Smith, 2008) with the luminosity-mass correlation by Salaris and Cassisi (2005). This gives the following relationship between main-sequence lifetime, τ MS , and stellar mass, M ⋆ , that is (Westby and Conselice, 2020),τ MS ≈ { 7.1 M ⋆ − 2.5 , 2.0 ≤ M ⋆ ≤ 20 M ⊙ 10 M ⋆ − 3 , 0.43 ≤ M ⋆ ≤ 2.0 M ⊙ 43 M ⋆ − 1.3 , M ⋆ ≤ 0.43 M ⊙ .(8)

By applying this relation, we obtain the black solid line in the upper panel of Figure 2, which illustrates the main-sequence lifetime of stars as a function of stellar mass. For comparison, the red line corresponds to a calculation performed with a similar relationship by Hansen and Kawaler (1994). The blue line additionally illustrates the “mean habitable lifetime for planets that may possess intelligent observers” by Waltham (2017) and was calculated with the power-law fit from the same study. Waltham (2017) based their mean habitable lifetime on a probabilistic combination of the Copernican and Anthropic principles and therefore assume planets to be inhabitable for intelligent life if these provide conditions similar to the Earth’s. Note, however, that (i) this lifetime is larger than the main-sequence lifetime for stars with masses below ∼0.6 M ⊙ , and (ii) that it is intended for intelligent observers and not for our definition of complex life (see Section 3.2); this lifetime should hence be taken with caution. The lower panel of Figure 2 additionally shows the fraction of stars that are still on the main sequence as a function of stellar age for the same relations.

OPEN IN VIEWER

FIG. 2.

Upper panel: Main-sequence lifetime of stars with masses between 0.1M ⊙ and 2.0 M ⊙ . The solid black line shows the lifetime as calculated through the approach by Westby and Conselice (2020) and as described in Equation 8. For comparison, the red line is calculated with a similar age–mass relationship as suggested by Hansen and Kawaler (1994), while the blue line shows the mean habitable lifetime for planets that may possess intelligent observers as proposed by Waltham (2017); note, however, that this estimate gives larger ages for low-mass stars than their respective main-sequence lifetime. Lower panel: The fraction of stars that are still on the main sequence as a function of age for the same three relations.

While we implement the relationship by Westby and Conselice (2020) as an upper limit of stellar age into our nominal case (the white area in Fig. 4 derives from this limit), we see in Section 5.2 that our results are insensitive to the chosen stellar main-sequence lifetime relation. This is because the upper age limit of stars, which allow for the existence of EHs, will always be set by the maximum bolometric luminosity that still permits the survival of complex life in the HZCL, a value that will always be lower than the actual main-sequence lifetime (but may intersect with the mean habitable lifetime by Waltham, 2017).

As a final step, we need to implement the actual SFH, that is, the evolution of the SFR of the galactic disk to obtain the distribution of presently existing main-sequence stars as a function of stellar mass and age. For this, we implement the SFH as reconstructed by Snaith et al. (2015) into our nominal case. These authors developed a chemical evolution model to reconstruct the SFH of the galactic disk over time from present-day chemical abundances of the Milky Way’s inner ( < 7 − 8 kpc) and outer disk. We chose this specific SFH by Snaith et al. (2015) over others, since their model also provides the metallicity evolution of the disk within one consistent framework. In addition, most of the other studies reconstruct a so-called cosmic SFH that is based on observations of SFRs in various galaxies throughout the Universe (van Dokkum et al., 2013; Madau and Dickinson, 2014; Stanway et al., 2018; Westby and Conselice, 2020). Recent studies by Ruiz-Lara et al. (2020) and Alzate et al. (2021), however, that reconstructed the SFH of the solar neighborhood based on the Gaia Data Release 2 (Gaia Collaboration et al., 2018), find an SFH that also shows bursts of star formation. The formation of the solar system took place during one of these major bursts that peaked at ∼ 5.5 Ga (Madau, 2023) and was likely triggered by an external perturbation (Mor et al., 2019) such as the first passage of the Sagittarius dwarf galaxy (Ruiz-Lara et al., 2020) or a nearby SN (Banerjee et al., 2016; Boss et al., 2008; Cameron and Truran, 1977).

The upper panel of Figure 3 shows the best-fit SFH for the inner (solid black line) and outer disks (dashed black line) by Snaith et al. (2015) that we use within our nominal case. For comparison, we also display the reconstructed SFH for the Milky Way by Naab and Ostriker (2006), the cosmic SFH (green line) by Madau and Dickinson (2014), and the SFH of a 2 kpc wide bubble around the Sun (gray line) by Ruiz-Lara et al. (2020). We evaluate how a change in SFH might affect the maximum number of EHs in Section 6.

OPEN IN VIEWER

FIG. 3.

Upper panel: The normalized best-fit star formation history (SFH) for inner (solid black line) and outer disk (dashed black line) by Snaith et al. (2015), the reconstructed SFH of the Milky Way by Naab and Ostriker (2006), the cosmic SFH (green line) by Madau and Dickinson (2014), and the SFH of a 2 kpc wide bubble around the Sun (gray line) by Ruiz-Lara et al. (2020). Lower panel: Inner and outer SFH by Snaith et al. (2015) renormalized (dotted lines) to only include stars that still reside on the main sequence.

The SFH in the upper panel of Figure 3, however, shows any star that ever existed in the galactic disk. Since we are only interested in those that are still on the main sequence, we need to renormalize the distribution. This can easily be done by excluding any star that already diverged from the main sequence via Equation 8. The lower panel of Figure 3 consequently shows the renormalized SFH for the inner and outer disk as implemented in our nominal case. The distribution, as expected, slightly shifts toward the present day.

As a side note, be aware that the SFH by Snaith et al. (2015) starts at an age of 14 Ga, even though the actual age of the Universe is currently estimated to be around 13.8 Gyr (Planck Collaboration et al., 2020) and the epoch of reionization with the emergence of the first stars might date some few 100 Myr later still (Robertson, 2021), with the oldest galaxies having emerged earlier than about 400 Myr after the Big Bang (Robertson et al., 2022).20 By using this specific SFH, we are stuck with this age determination, but this will not affect our results significantly. A later starting age would potentially decrease the maximum number of EHs further since the earliest stars would have had less time for their XUV flux to decrease (see Section 5.2).

By implementing IMF, SFH, and main-sequence lifetime into our model, we can already derive some interesting numbers even though we did not yet calculate N ⋆ itself (see also Table 6 for all input parameters).

From all the stars that ever formed within the disk, 91.15(+0.60/−2.04)%21 are yet residing on the main sequence. Within a stellar mass range of M ⋆ = 0.1 − 1.25 M ⊙ , even 97.39(+0.10/−0.0)% survived until the present day. For our nominal case with M ⋆ low = 0.08 M ⊙ , the IMF from Kroupa (2001), and the SFH from Snaith et al. (2015), only 0.45% of all presently existing stars have a mass that is >1.25 M ⊙ , while 12.78% are below 0.1 M ⊙ . Consequently, 86.77% of all stars in the disk reside between M ⋆ = 0.1 − 1.25 M ⊙ . For our maximum case with the IMF from Chabrier (2003) and the SFH from Naab and Ostriker (2006), none of the stars has a mass below 0.1 M ⊙ since M ⋆ low was chosen to be exactly 0.1 M ⊙ . In this study, 0.27% of all presently existing stars have a stellar mass that is >1.25 M ⊙ while 99.73% reside within M ⋆ = 0.1 − 1.25 M ⊙ . For our minimum case with M low = 0.07 M ⊙ , the IMF from Kroupa (2001), and the SFR from Snaith et al. (2015), 83.39% of all presently existing stars reside within M ⋆ = 0.1 − 1.25 M ⊙ while only 0.22% are above that mass range.

The average initial stellar mass calculated through the IMF is M ¯ ⋆ = 0.57(+0.08/−0.04) M ⊙ , while the average stellar mass for the present day is significantly smaller due to massive stars living much shorter, that is, M ¯ ⋆ = 0.29(+0.03/−0.01) M ⊙ . For comparison, Chabrier (2003) gives a present-day stellar number density in the disk of n ⋆ ∼ 0.134 pc^–3 and a mass density of ρ ⋆ ∼ 4.1 × 10 − 2 M e p c − 3 , which results in an average stellar mass of M ¯ ⋆ ∼ 0.31 M ⊙ . As initial values, these authors further give n ⋆ ∼ 0.145 pc^–3 and ρ ⋆ ∼ 8.2 × 10 − 2 M ⊙ p c − 3 equating to M ¯ e ∼ 0.57 M ⋆ 22.

Finally, only 46.86(+0.48/−3.33)% of the entire stellar mass that has ever been produced within the disk is yet residing on the main sequence, which is in stark contrast to the aforementioned stellar survival rate of 91.15(+0.60/−2.04)%. This illustrates a strong shift toward lower mass stars over time as can also be seen through a change in the relative stellar distribution within the different spectral classes. By considering the canonical IMF with M low = 0.08 M ⊙ , the respective fractions of M and K stars increase from their initial values of 72.74% and 14.28% to 79.80% and 15.67%, respectively. All the other spectral classes decrease in relative size, that is, the G and F stars from 3.76% and 2.96% to 3.38% and 0.84%, respectively. Higher spectral classes show an even stronger decrease due to their short lifetimes with present-day values of 0.23%, 0.07%, and <0.001% for A, B, and O stars, respectively.

The relative distribution of stars at the present day as a function of stellar birth age and mass can be seen in Figure 4 for our nominal case. The upper panel shows the relative fraction of stars within age and mass bins of δt = 50 Myr times δM = 0.01 M ⊙ and within M ⋆ = 0.1 − 1.25 M ⊙ . One can see that low-mass stars are dominating the distribution, and the ages of stronger SFRs from the implemented SFH by Snaith et al. (2015) can easily be identified. There are no stars in the white area at the upper right side of this figure since their main-sequence lifetime is too short to still exist at the present day. The lower panel depicts the fraction of stars born later than the respective age for all stars ever existing over the entire stellar mass range (dotted black line) and within M ⋆ = 0.1 – 1.25 M ⊙ (solid black line), respectively, and the same for today’s main-sequence stars (dotted and solid blue lines). There, it can be seen that only 23.0% of all main-sequence stars and 19.7% of those within M ⋆ = 0.1 – 1.25 M ⊙ are younger than the Sun.

OPEN IN VIEWER

FIG. 4.

Upper panel: The relative fraction of today’s main-sequence stars between 0.1 and 1.25 M ⊙ as a function of birth age and stellar mass within bins of 50 Myr and 0.01 M ⊙ for our nominal case. The varying densities at different ages are an effect of the implemented star formation rate (SFR) by Snaith et al. (2015) and depict higher and lower SFRs. The dashed orange lines indicate the birth age and mass of the Sun. Lower panel: The fraction of stars that are younger than the respective age, that is, for all stars ever existing over the entire stellar mass range (dotted black line), for stars between 0.1 M ⊙ and 1.25 M ⊙ (solid black line), and the same for today’s main-sequence stars (dotted and solid blue lines). Only 23% of stars between 0.1 M ⊙ and 1.25 M ⊙ are younger than the Sun. The green dashed line indicates an age of 6 Gyr, the approximate age geological activity will cease for cosmochemically Earth-like planets. A total of 71.5% of all stars are older than that age.

According to our nominal case, about 71.5% of all main-sequence stars are older than 6 Gyr. This age is of particular interest since data from coupled galactic chemical evolution and geophysical thermal models indicate that the window for geologically active, cosmochemically Earth-like planets is around 6 Gyr (Frank et al., 2014; Mojzsis, 2021). Older planets might therefore show strong restrictions in maintaining carbon–silicate and nitrogen cycles.

5.1.2. The metallicity requirement, α GHZ met

5.1.2.1. Calculating and implementing α GHZ met

Planet occurrence rates seem to be correlated with the metallicity of the host stars (Adibekyan, 2019; Zhu, 2019; Hansen et al., 2021). Based on observations, Fischer and Valenti (2005) were the first to present a probability relationship between a star’s metallicity, (Fe/H), and the occurrence rate of gas giant planets (GGPs) with orbital periods shorter than 4 years around FGK main-sequence stars, that is,P ( GGP ) = 0.03 × 10 2.0 [ Fe/H ] ,(12) where P (GGP) denotes a star’s probability of hosting such GGPs over the metallicity range − 0.5 < [ Fe/H ] < 0.5 . A correlation between giant planets and metallicity was later also confirmed by planet population synthesis models (Mordasini et al., 2012; Emsenhuber et al., 2021). In addition, Johnson and Li (2012) proposed a critical metallicity value needed for a planet to form that also depends on the orbital distance, r pl , from the star, that is,[ Fe/H ] crit ≃ − 1.5 + log ⁡ ( r pl / 1 A U ) .(13)

Based on this relation, Johnson and Li (2012) further suggested that the first Earth-like planets should have formed around stars with a metallicity of Z ⋆ ≳ 0.1 Z ⊙ .30 Planet population synthesis models (Emsenhuber et al., 2021) later found that the occurrence rate of Earth-like planets may peak at a metallicity around (Fe/H) = −0.2 since their formation is either limited by too few building blocks for lower metallicities or negatively affected for higher metallicities by more massive and potentially dynamic planets that either accrete or eject terrestrial planets. Earlier observations found an average metallicity for rocky exoplanet hosting stars of (Fe/H) = −0.02, that is, about solar (Buchhave et al., 2014). For low-mass stars, the same synthesis models obtained a strong positive dependence between metallicity and the occurrence rate of terrestrial HZ planets (Burn et al., 2021).

Zhu (2019) and Hansen et al. (2021) further investigated the relationship between planetary radius, R pl , and metallicity of the host star and found that the positive correlation between planet occurrence rate and metallicity still holds for small-sized planets below R pl < 4 R_⊕, where R_⊕ denotes the radius of the Earth. Based on the Kepler sample, all confirmed Kepler object of interest planets with such radii are found around stars with Z = Z ⋆ / Z ⊙ ≳ 0.25 , that is, with [ Fe/H ] ≥ − 0.6 (Adibekyan, 2019; Hansen et al., 2021). Jiang et al. (2021) further found a relationship between planet mass, M pl , and stellar metallicity; almost no planets with M pl ≲ 2 M_⊕ can be found for metallicities below [ Fe/H ] ∼ − 0.3 , that is, for Z ∼ 0.5 . However, predefining a definitive metallicity cutoff, Z min , for rocky exoplanets might at present be elusive. We therefore implement Z min = 0.3 for our nominal case, Z min = 0.1 for our maximum case, and Z min = 0.5 − 0.75 for our minimum case, respectively.

To implement the importance of metallicity (i.e., α GHZ met ) into our model, we have to consider

the spatial metallicity distribution within the present galactic disk (see Fig. 7, upper panel), and

the temporal metallicity evolution of the galactic disk over time (see Fig. 7, lower panel).

If we did not account for the evolution of galactic metallicity, we would overestimate G- and F-type stars with lower metallicities but at the same time underestimate M and K dwarfs with lower metallicities, since old, low-metallicity stars are predominantly low-mass stars.

To take account of the present-day metallicity, we closely follow Westby and Conselice (2020) and implement the galactic MDFs suggested by Hayden et al. (2015). These authors provide skewed Gaussian metallicity distributions, including mean values, standard deviation, and the respective skewness, for different galactocentric distances within the disk. Similar to Westby and Conselice (2020), we use these values to calculate the galactic probability distribution of stars above a certain metallicity threshold. Since Hayden et al. (2015) only provide data between 3 kpc and 15 kpc, we make use of their fitting relationship for the inner and outer skirts of the disk, that is, [ Fe/H ] = − 0.053 × r + 0.428 . For z = 2.0 − 2.5 kpc, we assume the same mean metallicity value as given by Hayden et al. (2015) for z = 1.0 − 2.0 kpc, which may slightly overestimate the number of stars above our metallicity thresholds. Similarly, skewness and standard deviations are assumed to stay constant outward into the disk. The upper panel of Figure 7 shows the mean values for the metallicity distribution according to Hayden et al. (2015) with the extensions discussed in this paragraph. Within r = 3 − 15 kpc and z = 0 − 2 kpc, the range for which these authors provide observational data, the highest mean value for the metallicity of [ Fe/H ] = + 0.11 was measured at a galactocentric distance and height of 5 < r < 7 kpc and 0 < z < 0.5 kpc, respectively, while the lowest measured value of [ Fe/H ] = − 0.42 was found surprisingly close to the center at 3 < r < 5 kpc but for 1 < z < 2 kpc.

OPEN IN VIEWER

FIG. 7.

Upper panel: The metallicity (Fe/H) distribution of the present-day galactic disk according to Hayden et al. (2015). Closer than r = 3 kpc and beyond r = 15 kpc the value for (Fe/H) was extrapolated with the fitting function provided by Hayden et al. (2015); for z > 2 kpc it was further assumed to be the same as for z < 2 kpc. Lower panel: Evolution of (Fe/H) for the inner (red; model with dilution; see text) and outer disk (blue) according to the best-fit model by Snaith et al. (2015).

As we further account for the metallicity evolution of the Milky Way, we implement the best-fit (Fe/H)-age tracks provided by Snaith et al. (2015) for the inner and outer disk, which correlate with the best-fit SFH from their chemical evolution model. These tracks can be seen in the lower panel of Figure 7 where “with dilution” denotes their best-fit model for the outer disk. For this, Snaith et al. (2015) assumed a galactic accretion event at 10 Ga that diluted the in situ gas of the outer disk with primordial gas to match the age-metallicity and age-α observations by Haywood et al. (2013). It is worth noting that a similar dilution effect was later also considered by Spitoni et al. (2019, 2020, 2021) and Lian et al. (2020b, 2020a) within their chemical evolution models for reproducing the abundance ratio of (α/Fe)31 versus (Fe/H) in Apache Point Observatory Galactic Evolution Experiment32 stars in the outer disk.

To combine the evolutionary tracks from Snaith et al. (2015) with the present-day distribution from Hayden et al. (2015), we renormalize the (Fe/H)-age tracks for inner and outer disk with the present-day values of the different Milky Way regions so that the evolving metallicity of the best-fit model by Snaith et al. (2015) will reach the certain regional (Fe/H)-values as observed by Hayden et al. (2015). The average metallicity of all stars within a certain region will hence be renormalized to meet the present-day metallicity value of this respective region. If we generally average over the (Fe/H)-values provided in Hayden et al. (2015) for the outer disk region between 7 and 15 kpc (for a disk height 0 < z < 0.5 kpc), we receive (Fe/H) = −0.19, which is close to (Fe/H) approximately −0.20, the value Snaith et al. (2015) receive for their outer disk for the present day (see their Fig. 9d). For the inner disk, Snaith et al. (2015) find (Fe/H) approximately +0.35 which is slightly higher than (Fe/H) = +0.23, the average value in Hayden et al. (2015) for 3 < r < 7 kpc and 0 < z < 0.5 kpc. Note, however, that the inner disk in Snaith et al. (2015) covers the entire disk with radii r < 7 − 8 kpc, which may explain such a higher value. To calculate the Gaussian distribution, we assume that mean, standard, and skewness remained constant over the entire age of the Milky Way.

5.1.2.2. The effect of α GHZ met

The upper panel of Figure 8 shows the probability of any presently existing star to be above the threshold value of our nominal case, that is, Z min = 0.3 , as a function of galactocentric distance r and height z and based on today’s galactic (Fe/H)-values observed by Hayden et al. (2015). As described above, each region33 with a different present-day (Fe/H) value will have a slightly different (Fe/H) evolution over time, as we require that the average metallicity of all stars presently existing within such a region has to reach the present-day value of (Fe/H) finally. The middle panel of Figure 8 therefore illustrates the probability of being above such threshold for each time step averaged over all galactic regions for the inner (red line) and outer disk (blue line). One can see that the galactic accretion event at 10 Ga significantly reduced the probability of a star being above the metallicity threshold. The light red and blue lines further show those galactic regions in the inner and outer disks that presently have the minimum and maximum probability values, respectively.

OPEN IN VIEWER

FIG. 8.

Upper panel: The present-day probability for a certain star to be above a metallicity threshold of Z min = 0.3 (nominal case) as a function of galactocentric distance r and height z. This probability is based on the metallicity distribution functions by Hayden et al. (2015). Middle panel: The probability of being above the same threshold ( Z min = 0.3 ) for each time step for the inner (nontransparent red line) and outer disk (nontransparent blue line). The probability is averaged over each region (i.e., spatial bin within our model) that holds a different present-day value of (Fe/H) within the inner and outer disk, respectively. The probabilities to be above Z min = 0.3 of those regions with the lowest and highest (Fe/H)-values in the inner and outer disk are additionally visualized by the transparent red and blue lines, respectively, to illustrate the entire probability spread in the disk. Here, the probability is based on the evolutionary (Fe/H)-tracks from the best-fit model by Snaith et al. (2015) but renormalized to the present-day values of each region. Lower panel: The mass-age distribution of stars within the GHZ that have a metallicity of Z > 0.3 (within a bin size of 0.1 M ⊙ times 50 Myr). As can be seen, most of the oldest stars do not fulfill the metallicity requirement, which the middle panel can easily explain.

The lower panel of Figure 8 illustrates the age distribution of remaining stars that fulfill our nominal case metallicity threshold of Z min = 0.3 . There is a strong increase in stars above the threshold immediately after the galaxy’s birth as many of the oldest stars have a too low metallicity to allow the formation of rocky exoplanets. However, because the (Fe/H)-tracks are steeply rising (as seen in this Fig.’s middle panel), this effect is mostly restricted to early ages. At the Sun’s birth age, for instance, the mean metallicity of stars born at the same age and region measures [ Fe/H ] ∼ 0.12 , a value that is above the solar value.

In total, 69.3% of all remaining stars meet the metallicity threshold of Z min = 0.3 , that is, 1.65 × 10 10 stars. For Z min = 0.1 , this value rises toward 2.42 × 10 10 stars, while it decreases toward 1.30 × 10 10 stars for Z min = 0.5 . If we increase the metallicity threshold to Z min = 0.75 , only 0.94 × 10 10 stellar candidates remain. By comparing the latter result with GHZ studies that include the effect of metallicity, this value closely resembles the range of stars we obtain by implementing the probability distributions of Spitoni et al. (2014) and Spitoni et al. (2017) for which we find 0.70 × 10 10 to 1.16 × 10 10 stars. By implementing the one by Vukotić et al. (2016), we obtain ∼ 2.0 × 10 10 stars, a value that is relatively close to 2.42 × 10 10 stars for our case with Z min = 0.1 . It is therefore reasonable to take Z min = 0.1 as our maximum and Z min = 0.75 as our minimum case.34 If we also vary SFH, IMF, etc., we obtain a total range of 1.65 ( + 0.77 / − 0.72 ) × 10 10 stars.

Finally, 76.69(+1.00/−0.45)% of the remaining stars are low-mass M dwarfs, which is a slightly lower value than for the entirety of GHZ stars, for which it is 77.27(+0.77/−0.0)%. This small shift relates to the fact that the metallicity of the ISM has been increasing since the galaxy’s birth, thereby implying that M dwarfs are, on average, below the threshold more often than heavier stars. This can be seen in Figure 9, where the upper panel illustrates the stellar fraction above the respective metallicity threshold Z min . This fraction remains constant for stars with a main-sequence lifetime larger than the present age of the galaxy. As soon as the main-sequence lifetime becomes shorter, however, the stellar fraction starts to increase since stars born close to the birth age of the Milky Way already left the main sequence. This is also exemplified in this figure’s lower panel, which shows the number of stars above the respective Z min , again as a function of stellar mass. For larger Z min , more M and K dwarfs are below the threshold and the relative fraction of G and F stars is consequently rising. A total of 17.83(+0.10/−1.54)% of the remaining stellar sample are hence K dwarfs, while 4.43(+0.34/−0.17)% and 1.15(+0.34/−0.03)% are G and F stars, respectively.

OPEN IN VIEWER

FIG. 9.

Upper panel: The fraction of stars above the different metallicity thresholds Z min as a function of stellar mass. Lower panel: The number of stars above Z min , again as a function of stellar mass.

As a final note related to metallicity, we point out that the occurrence rate of rocky exoplanets will not only depend on a lower metallicity cutoff. Their formation probability may also decrease for high metallicities. This seems to be likely due to (i) a rise in the occurrence rates of hot Jupiters (Buchhave et al., 2018; Osborn and Bayliss, 2019) and warm, dynamically active giant planets (Schlecker et al., 2021) with increasing metallicity, and (ii) an increase in rocky building blocks that may lead to the growth of more massive planets with hydrogen-dominated atmospheres. Through assuming (i) and the related destruction of rocky exoplanets by migrating hot Jupiters, Prantzos (2008) and subsequently Carigi et al. (2013), used Equation 12 by Fischer and Valenti (2005) to calculate the probability of a star hosting Earth-like planets but not hot Jupiters, thereby taking into account a potential lower probability of rocky exoplanets around high-metallicity stars. However, Spitoni et al. (2014) already pointed out that this assumption might be questionable since the relationship by Fischer and Valenti (2005) covers any GGPs with orbital periods shorter than 4 years.35 Since GGPs on wider orbits (as is the case in the solar system) will not necessarily destroy or disturb Earth-like planets in circumstellar HZs, one may expect that applying such relationship could slightly underestimate the occurrence rate of rocky exoplanets around high-metallicity stars (see also Fig. 8 in Carigi et al., 2013), whereas the opposite may hold for neglecting any such relationship (as in our model). This reasoning, however, is only based on argument (i), while the relevance of argument (ii) cannot be properly assessed with our present scientific knowledge.

5.2.1. The lower stellar limit, α at XUV

5.2.1.1. Implementing α at XUV

A star’s XUV irradiation is absorbed in the upper atmosphere of a planet, which leads to the heating and expansion of its thermosphere (see Fig. 35 in the Appendix). As we detail in the appended Section 1.2, this effect can be specifically important for Earth-like atmospheres, for which XUV surface fluxes only 5 to 6 times as high as received by today’s Earth will be enough for the atmosphere to expand adiabatically and to erode into space within a few Myr (Johnstone et al., 2021b; Tian et al., 2008a). Since CO₂ serves as an infrared coolant (Johnstone et al., 2021b; Kulikov et al., 2007; Lichtenegger et al., 2010), it can prolong the thermal stability of Earth-like atmospheres and we derive certain stability thresholds for nitrogen-dominated atmospheres with maximum mixing ratios of x CO 2 , max = 1% and x CO 2 , max = 10% CO₂ in the appended Section 1.2. This important phenomenon implies that stars with certain XUV surface fluxes in the HZCL cannot host planets with stable N₂-O₂-dominated atmospheres.37 The lower stellar limit, α at ll of A at ( ∏ i = 1 n α at i ) , is therefore of crucial importance for the distribution of EHs in the galaxy.

For implementing α at ll , which we simply implement as XUV limit, α at XUV , we calculate the X-ray and/or XUV surface flux, F X and F XUV , in the middle of the HZCL for a range of stellar masses. For this, we make use of the Python package Mors,38 which solves the stellar evolution model developed by Johnstone et al. (2021b). The package implements a grid of stellar isochrones by Spada et al. (2013) as input for the basic physical parameters and it can be used to simulate derived stellar parameters such as XUV luminosities as functions of stellar mass and age for any percentile of the stellar rotational distribution.39 Noteworthy, Kubyshkina and Fossati (2022) were able to reproduce the mass-radius distribution of presently known exoplanets as a function of stellar activity and evolution by implementing Mors into their model. Since Mors can simulate any star between M ⋆ = 0.1 − 1.25 M ⊙ , this also predefines the stellar mass range we consider within our model. Although this approach excludes any low-mass M dwarfs below M _å = 0.1 M ⊙ such as Trappist-1, this lower limit does not significantly affect our derived results since, among other reasons (see Section 7.4), Earth-like atmospheres are unlikely to be thermally stable around these low-mass dwarfs, as derived in the following sections.

Through this approach, we exclude any star for which F X and/or F XUV will be too high to permit the long-term stability of an N₂-O₂-dominated atmosphere. This statement, however, deserves some further explanation.

We focus on the HZCL since this concept provides important boundary conditions for our model. Schwieterman et al. (2019a), and more recently Ramirez (2020), define the HZCL (or CLHZ as abbreviated by Ramirez, 2020) as the zone around a star where (i) liquid water can exist at a planet’s surface and (ii) the CO₂ partial pressure is below the toxicity limit for complex life (see Section 3.2.5). Both studies find a partial pressure of pCO 2 ∼ 0.1 bar to be the upper toxicity limit for advanced metazoans, that is, all animals except sponges.40 This pressure limit refines the outer HZ boundary since the farther away a planet might be from its host star, the more greenhouse gases will be needed to counterbalance the decrease in the stellar surface flux. A maximum amount of pCO₂ therefore restricts such boundary and defines the HZCL.

Schwieterman et al. (2019a) provide HZCL boundaries for three different values of the CO₂ partial surface pressure, that is, pCO 2 = { 0.01 , 0.1 , 1.0 } bar, and show that for a solar-like star, the HZCL is only 21%, 32%, and 50% as wide as the conventional HZ, respectively. Ramirez (2020) refines the HZCL for pCO₂ = 0.1 bar by using an advanced energy balance model that also includes the effect of clouds and finds that it would be ∼35% wider than initially found by Schwieterman et al. (2019a). Trace species such as N₂O and CO₂, pressure broadening by N₂, planetary and atmospheric mass, continental distribution, surface gravity, the surface water content, or planetary rotation rate, however, can further alter HZ or HZCL boundaries (Kodama et al., 2019; Kopparapu et al., 2014; Schwieterman et al., 2019a; Vladilo et al., 2013).

Since we are interested in the maximum number of EHs, that is, HZCL planets that can host N₂-O₂-dominated atmospheres with minor amounts of CO₂, we can combine the results obtained by Johnstone et al. (2021b) on the stability of N₂-dominated atmospheres with the concept of the HZCL. For this, we assume any planet to be in the middle of the HZCL even though the actual positions of HZCL planets will certainly vary over its entire orbital range. However, planets inside the mean HZCL distance, d 〈 HZCL 〉 , will receive a correspondingly higher amount of F XUV . In contrast, planets farther out will receive a correspondingly lower one. By assuming that all the HZCL planets are evenly distributed over the entire HZCL, we can simplify the problem to the mean value, d 〈 HZCL 〉 . A smaller HZCL will further result in d 〈 HZCL 〉 being closer to the star.

Similarly, we assume that all planets have a mass of M _pl = 1.0 M_⊕ since this is the mass for which reliable simulations are available. Planets with higher masses than the Earth will tend to grow faster to reach such a high mass, and the more massive a planet, the likelier it will become that any accreted primordial atmosphere cannot be lost afterward. The so-called Fulton Gap (Fulton et al., 2017; Fulton and Petigura, 2018, see also Section A.2) further supports the hypothesis that planets of a certain mass/radius will tend to host a primordial atmosphere while planets with a radius of R pl ≲ 1.5 R_⊕ (Rogers, 2015), or even as low as R pl ≲ 1.23 R_⊕ (Chen and Kipping, 2017), likely end up being rocky exoplanets,41 at least for orbital periods with P pl < 100 days (Fulton et al., 2017). Such a radius would give an upper limit for the mass of a rocky exoplanet of M pl ≲ 4 M_⊕. However, for planets orbiting farther out from their host star, this value might become significantly lower since photoevaporation will be substantially weaker in the HZCL of solar-like stars. For a planet that orbits a G-type star at 1 AU, an accreted H₂-dominated atmosphere might not be eroded subsequently, if the planet grows to a mass of only M pl ∼ 1 M_⊕ within the protoplanetary gas disk (Erkaev et al., 2023). This mass, however, will certainly be higher for lower mass stars, and there are indications that planets around more massive stars indeed tend to have larger primordial atmospheres (Lozovsky et al., 2021).

Another upper mass limit can be derived from the degassing of volatiles from a planet’s interior and, potentially, its tectonic regime. While Valencia et al. (2007) and van Heck and Tackley (2011) find that plate tectonics might be equally or more frequent on super-Earths, several other studies conclude that the opposite might be true (Miyagoshi et al., 2015; Noack and Breuer, 2014; O’Neill and Lenardic, 2007). In addition, Noack et al. (2017) and Dorn et al. (2018) for stagnant-lid, and Kruijver et al. (2021) for plate tectonic planets, found that the degassing of volatiles will start to be significantly reduced for planets with M pl ≳ 3 − 4 M_⊕. This provides another indication that the upper planetary mass limit may be found below M pl ∼ 4 M_⊕. A further constraint on the upper mass limit for habitable planets can be derived from the lifetime of a primordial continent. Dohm et al. (2018) argue that a primordial anorthositic continent on super-Earths will quickly be transported into the interior due to intense mantle convection, thereby strongly limiting the supply of nutrients needed for the emergence of life. If so, the origin of life would have to happen earlier on heavier planets than on lighter, Earth-like, ones.

On the other side of the mass scale, small-mass planets will not be able to hold on to an N₂-dominated atmosphere. A Mars-mass planet around a G-type star might already be too small to assure the long-term stability of such an envelope and even more massive planets might struggle to keep such an atmosphere, particularly for K and M dwarfs. In fact, nonthermal atmospheric escape tends to reach maximum loss rates for a planetary radius of R pl ∼ 0.7 R_⊕ (Chin et al., 2024),42 indicating that the lower radius of a habitable planet may approach R pl ≳ 0.7 R_⊕. Small-mass planets will also need longer to provide stable conditions for an N₂-O₂-dominated atmosphere against escape into space but, at the same time, will become geologically inactive within a shorter duration. At some age, these limits will start to overlap, and these will overlap faster for small-mass planets than for higher mass ones.

These arguments exemplify that for both sides of the mass distribution, planets will become inhospitable to complex life for certain mass limits. However, these upper and lower limits are not (yet) strictly definable, and their mean value might diverge from M pl = 1.0 M_⊕. However, because the distribution of planetary masses will likely have a peak value between log M pl / M ⊕ = 0.6 − 1.0 (Malhotra, 2015), assuming a mean value of 1 M_⊕ might potentially be an overestimate, as long as the minimum mass is not close to or at 1 M_⊕.

As a first step for implementing the XUV limit, α at XUV , however, we need to calculate the HZ boundaries. For this, Kopparapu et al. (2013) give the following relationship between the stellar surface flux, S eff , and the stellar effective temperature, T eff , that is,S eff = S eff , ⊙ + a T ⋆ + b T ⋆ 2 + c T ⋆ 3 + d T ⋆ 4 ,(14) where T ⋆ = T eff − 5780 K. The coefficients are dependent on the chosen inner and outer boundary, which typically resemble limits such as Water Loss, Recent Venus, or Runaway Greenhouse for the inner, and Moist Greenhouse, Maximum Greenhouse, and Recent Mars for the outer boundary (see Kasting et al., 1993; Kopparapu et al., 2013, 2014, for a discussion). The parameter S eff , ⊙ further gives the solar effective surface flux at the respective HZ boundaries in units of the solar surface flux, S eff , ⊕ , that reaches the top of today’s Earth’s atmosphere.

We vary between Moist Greenhouse and Runaway Greenhouse limits as the inner boundary, which coincide with S eff , ⊙ ∼ 1.014 S eff , ⊕ and S eff , ⊙ ∼ 1.0512 S eff , ⊕ , respectively, in Kopparapu et al. (2013), and with S eff , ⊙ ∼ 1.107 S eff , ⊕ for the Runaway Greenhouse in Kopparapu et al. (2014). At the Moist Greenhouse limit, the planetary atmospheric temperature increases to such an extent ( T pl ≳ 330 K) that the stratosphere becomes dominated by H₂O and an Earth-ocean can be lost on Gyr timescales (Kasting et al., 2015; Kasting et al., 1993; Kopparapu et al., 2013; Wolf et al., 2017). At the Moist Greenhouse limit, on the contrary, the Sun/star-induced heating reaches a certain threshold at which the thermal equilibrium can only be maintained by the evaporation of the entire surface water (Goldblatt et al., 2013; Leconte et al., 2013b; Nakajima et al., 1992). These boundaries were calculated by Kopparapu et al. (2013) to be at 0.99 AU and 0.97 AU, respectively, for the present-day Sun.

We substitute the outer boundary of the HZ with the HZCL definitions by Schwieterman et al. (2019a) and Ramirez (2020). While Schwieterman et al. (2019a) give the corresponding values for a, b, c, and d, we directly take S eff for a 1 bar atmosphere from Figure 2 in Ramirez (2020). The parameter T eff we further calculate with Mors for stellar masses between M ⋆ = 0.1 − 1.25 M ⊙ and for the entire lifetime of the respective stars (with the present age of the Universe, that is, the Hubble time of t 0 ∼ 13.8 Gyr, as maximum age). The corresponding HZ distances can then be calculated via (Kopparapu et al., 2013, 2014) the following:d HZ ( CL ) = ( L bol / L bol , ⊙ S eff ) 0.5 ,(15) where L bol , ⊙ is the solar and L bol the stellar bolometric luminosity, of which the latter will also be simulated through Mors, the Python package provided by Johnstone et al. (2021a). The mean HZCL distance, d 〈 HZCL 〉 , for which we calculate our threshold, is then simply defined as the average between the inner and outer HZCL boundaries.

For our nominal case, we take the outer HZCL boundary by Ramirez (2020). By additionally taking the Runaway Greenhouse limit by Kopparapu et al. (2014) as the inner boundary, this results in a mean HZCL distance of d 〈 HZCL 〉 = 1.17 AU for a solar-like star at 5 Gyr. Figure 10 further shows several different HZ boundaries for an age of 5 Gyr calculated with Mors and Equations 14 and 15. The half-transparent lines in the figure’s background illustrate the same boundaries but for an age of 100 Myr. Note that for 5 Gyr the Earth position (dashed black lines) is already slightly outside the HZCL. This is because the bolometric luminosity of the stellar evolutionary track that closest resembles the Sun in the model of Spada et al. (2013) is slightly higher at 5 Gyr than for the real Sun (i.e., L bol / L bol , ⊙ = 1.077 ).

OPEN IN VIEWER

FIG. 10.

The inner habitable zone (HZ) boundary (thick gray line) of the HZ according to the Runaway Greenhouse threshold and several different outer boundaries, that is, the HZ of complex life (HZCL) for pCO 2 = 0.1 bar (green line) and pCO 2 = 0.01 bar (green dashed line), both according to Schwieterman et al. (2019a), for pCO 2 = 0.1 bar (blue line) according to Ramirez (2020), and for the Maximum Greenhouse threshold (black line). All these boundaries are calculated through Mors and are for an age of 5 Gyr. The respective light colors in the background illustrate the same HZ boundaries but for 100 Myr. The dashed black lines illustrate the position of the Earth.

To be consistent with the outer HZCL boundary in Ramirez (2020), we assume a maximum of pCO 2 = 0.1 bar for the atmosphere. As can be seen in Figure 35 of appended Section 1.2, this threshold coincides with a maximum value for the X-ray surface flux at d 〈 HZCL 〉 of F X ∼ 5 erg s^–1 cm^–2 or F X ∼ 7.7 F X , ⊕ ,43 respectively (vertical solid black line). For such a flux, the exobase of an N ₂-O₂-dominated atmosphere with a partial pressure of pCO 2 = 0.1 bar (solid violet line) is already expanded to such an extent that a complete erosion of the atmosphere will take place within less than ∼10 Myr Johnstone et al. (2021b). In this study, it has to be noted that the simulations by Johnstone et al. (2021b) were performed with reconstructed synthetic solar spectra. As outlined in more detail in appended Section 1.1, however, the Sun seems to be particularly weak in the X-ray wavelength part of the spectrum but slightly less weak in its EUV part when compared with other stars (Johnstone et al., 2021a). Only considering F X as a threshold might, therefore, underestimate the number of stars that can presently hold an N₂-O₂-dominated atmosphere since the ratio of F X / F EUV is on average lower for the entire rotational distribution of stellar masses between 0.1 and 1.25 M ⊙ than for the Sun. To counter this effect, we take into account the whole XUV wavelength range (i.e., F XUV = F X + F EUV ) for our nominal case, which results in a threshold value for the XUV surface flux of F XUV , max = 35 erg s^–1 cm^–2 for a wavelength range of 0.5–92 nm and based on a ratio of F EUV = 6 F X (Tu et al., 2015; Johnstone et al., 2021a). In the lower panel of Figure 34, one can see that such a threshold value coincides well with the XUV flux of a slow rotator for the same age of ∼600 Myr. One can, therefore, expect that this threshold value might even be slightly “optimistic” because F X / F EUV is smaller for the Sun than for the average stellar distribution.

For our minimum case, we, consequently, take into account only the X-ray part of the spectrum with F X , max = 5 erg s^–1 cm^–2. For the maximum case, on the contrary, we assume F X , max = 8 erg s^–1 cm^–2 (vertical dashed line in Fig. 35 and the entire XUV spectrum with F XUV , max = 56 erg s^–1 cm^–2. As can be seen in Fig. 35, not even an N₂-O₂ atmosphere that contains 25% CO₂ would be stable for these conditions according to the results of Johnstone et al. (2021a). However, such an optimistic case might cover any potential parameters that might generally enhance the stability of an N₂-dominated atmosphere such as unknown cooling agents (additional atomic line cooling by oxygen as in Nakajima et al., 2022), an average mass for habitable rocky exoplanets that is significantly higher than 1.0M_⊕, and/or the unlikely case that the actual short-wavelength activity of the stellar distribution is significantly weaker than found by Johnstone et al. (2021a).

5.2.1.2. The effect of α at XUV

Figure 11 shows the X-ray (upper panel) and XUV surface flux (lower panel) at d 〈 HZCL 〉 for stars at the 50th percentile of the stellar rotational distribution as a function of stellar mass and birth age, calculated through the stellar evolution model Mors from Johnstone et al. (2021a). The orange contour lines illustrate the threshold values of our minimum (upper panel) and nominal cases (lower panel) for atmospheres with x CO 2 , max = 10%, that is, F X , max = 5 erg s^–1 cm^–2 and F XUV , max = 35 erg s^–1 cm^–2, respectively. All stars that are placed above the respective contour lines have a surface flux that is below the associated threshold value, implying that under such assumption, an N₂-O₂-dominated atmosphere with 10% CO₂ is stable. In this study, it can be seen that a higher number of stars are fulfilling the F XUV , max -levels than the related F X , max -levels, which stems from the finding by Johnstone et al. (2021a) that the Sun’s F X / F EUV ratio is lower than the average stellar ratio. The violet contour line in the lower panel further illustrates our maximum case for x CO 2 , max = 10 % with F XUV , max = 56 erg s^–1 cm^–2 and its respective value for the X-ray surface flux of F X = 8 erg s^–1 cm^–2 can be seen in the upper panel. The gray lines show the threshold values of our minimum case for x CO 2 , max = 1 % . All threshold values for α at XUV are listed in Table 7 of Section 8.

OPEN IN VIEWER

FIG. 11.

Upper panel: The stellar X-ray surface flux, F X , in the middle of the HZCL for moderately rotating stars (i.e., the 50th percentile of the stellar rotational distribution) with stellar masses between M ⋆ = 0.1 − 1.25 M ⊙ as a function of stellar mass and age, calculated with Mors (Johnstone et al., 2021a) and based on a grid of stellar isochrones by Spada et al. (2013). The plot shows three colored contour lines with different threshold values for the X-ray surface flux, that is, F X , max = { 2.5 , 5.0 , 8.0 } erg s − 1 cm − 2 . Lower panel: The same but for the X-ray and extreme ultraviolet (XUV) surface flux with colored contour lines at F XUV , max = { 17.5 , 35.0 , 56.0 } erg s − 1 cm − 2 .

It has to be noted, however, that Figure 11 only illustrates the 50th percentile of the stellar distribution, that is, the so-called moderate rotators. Stars with the same mass but at a lower percentile of the rotational distribution will rotate slower (such as the so-called slow rotator at the 5th percentile) and therefore reach any surface flux threshold earlier than the 50th percentile. Conversely, any star that rotates faster will fall below such value later. To account for this effect, we implement the entire stellar rotational distribution, as it can be calculated with Mors, and distribute all stars evenly over all percentiles.

Figure 12 illustrates the age at which each percentile (from 0 to 100th) of each stellar mass (with Δ M ⋆ = 0.01 M ⊙ ) will decline below the threshold values of both our nominal and minimum cases for x CO 2 , max = 10 % with F XUV , max = 35.0 erg s − 1 cm − 2 (lower panel) and F X , max = 5.0 erg s − 1 cm − 2 (upper panel), respectively. In this study, two features are particularly noteworthy. First, earlier percentiles, that is, slower rotating and therefore less active stars of the same stellar mass, fall below the given threshold values earlier. However, this effect is more pronounced for higher mass stars (see also Section A.1 and Fig. 34). Second, not all stars will meet such criteria during the Hubble time, t ₀, but F XUV , max will be met slightly faster and by a higher number of stellar masses than F X , max . For the given thresholds in Figure 12, stars with masses of M ⋆ ≥ 0.44 M ⊙ and M ⋆ ≥ 0.51 M ⊙ will fall below F XUV , max and F X , max , respectively, within the first 12 Gyr of their lifetimes. This already illustrates that most M dwarfs will never meet the criteria for our nominal and minimum cases since their short wavelength radiation declines far too slowly to provide habitable conditions within the current lifetime of the Universe.

OPEN IN VIEWER

FIG. 12.

The stellar ages at which the X-ray (upper panel) and XUV surface flux (lower panel) at d 〈 HZCL 〉 will fall below threshold values of (upper panel) and (lower panel), respectively, for each percentile (from 0 to 100th) of the stellar rotational distribution and stellar masses of = 0.1–1.25 M ⊙ . The greenish to reddish lines indicate the different stellar masses (from 1.25 M ⊙ ) with masses declining toward the right side of the plot with a step size of M ⊙ . For F X , max = 5.0 erg s − 1 cm − 2 , the threshold value can only be met within 12 Gyr for stars that have a mass of M ⋆ ≥ 0.51 M ⊙ while can be met already for M ⋆ ≥ 0.44 M ⊙ .

With the Runaway Greenhouse limit from Kopparapu et al. (2014) and the pCO₂ = 0.1 bar limit from Ramirez (2020) as the inner and outer HZCL boundaries, respectively, and by taking account of the entire stellar rotational distribution, we find that 2.44 × 10 9 stars out of the initial 1.64 × 10 10 GHZ stars fulfill the XUV requirement, α at XUV , for our nominal case with x CO 2 , max = 10 % and F XUV , max = 35.0 erg s − 1 cm − 2 . This is only ∼ 15.4 % and 3.1% of all remaining stars and all disk stars, respectively. By only changing the X-ray or XUV surface flux to F X , max = 5.0 erg s − 1 cm − 2 (minimum case) and F XUV , max = 56.0 erg s − 1 cm − 2 (maximum case), we obtain 1.65 × 10 9 and 4.54 × 10 9 stars, respectively. By varying all input parameters, including the inner HZCL boundary, which we change to the Runaway Greenhouse limit by Kopparapu et al. (2013) for the minimum case, we obtain an entire range of 2.44 ( + 3.93 / − 1.79 ) × 10 9 stars. For our minimum case, this indicates that only 0.79% of all disk stars (i.e., 0.65 × 10 9 stars) are still fulfilling all implemented necessary requirements needed to host an EH.

As a side note, if we change the outer HZCL boundary of our nominal case from Ramirez (2020) to Schwieterman et al. (2019a), the number of suitable stars decreases from 2.44 × 10 9 to 2.20 × 10 9 stars due to the HZCL being smaller and d 〈 HZCL 〉 being closer to the star. If we take the Moist Greenhouse limit from Kopparapu et al. (2013) instead of the Runaway Greenhouse limit from Kopparapu et al. (2014), but keep everything else within our nominal case the same, we obtain a slightly larger fraction of suitable stars, that is, 2.52 × 10 9 .

Figure 13 shows the number of remaining GHZ stars that are above the lower limit, α at XUV , for our nominal (upper), minimum (middle), and maximum (lower panel) cases. Even though these cases do not only have varying values for α at XUV and also minimum and maximum values for all other input parameters, the effect of the lower limit can easily be seen. The lower the threshold values for F XUV , max or F X , max are, the higher the stellar mass must be to still meet this requirement. While for our nominal case, all stars below M ⋆ = 0.41 M ⊙ fail to meet this necessary criterion, this threshold changes toward M ⋆ = 0.22 M ⊙ and M ⋆ = 0.50 M ⊙ for the maximum and minimum case, respectively. Whether M dwarfs are still part of the remaining stellar distribution, therefore, strongly depends on the lower limit. Stellar masses below M ⋆ ∼ 0.20 M ⊙ , however, will likely never meet this requirement, at least as long as their XUV flux cannot be considered exceptionally weak. This also exemplifies why a priori excluding stars with stellar masses below M ⋆ = 0.1 M ⊙ can at most have negligible effects on our results.

OPEN IN VIEWER

FIG. 13.

The number of remaining GHZ stars that are above the lower limit, α at XUV , as a function of age and stellar mass within a bin of Δ M ⋆ = 0.01 M ⊙ and Δ t = 50 Myr. The upper panel shows our nominal case with F XUV , max = 35.0 erg s − 1 cm − 2 , and the outer HZCL boundary from Ramirez (2020) with our standard values for the initial mass function (IMF), SFR, and GHZ. The middle and lower panels illustrate our minimum and maximum cases with F X , max = 5.0 erg s − 1 cm − 2 and F XUV , max = 56.0 erg s − 1 cm − 2 , respectively, combined with the minimum and maximum values of all other input parameters.

The lower limit therefore strongly affects the relative fractions of the different spectral classes, much more so than the metallicity threshold. Only 0–22.36% of the remaining stars will be M dwarfs, while K dwarfs suddenly hold the majority with 53.38–68.66%. The relative fraction of G and F stars also rise significantly, that is, from 4.26–4.76% and 1.12–1.49% to 18.09–35.54% and 4.36–11.08%, respectively.

Figure 14 illustrates this distribution. In this study, the upper panel depicts our nominal case and shows the number of stars by spectral class that still meet all requirements as a function of birth age. At maximum, only old M dwarfs remain part of the sample (for the minimum case, no M dwarfs remain, not even old ones), while the distribution shifts toward younger ages for earlier spectral classes. This is particularly interesting because any N₂-dominated atmosphere must be outgassed at an age when the star already allows for its thermal stability. A geologically active planet must therefore be a prerequisite for allowing the buildup of a stable atmosphere. As already briefly discussed in Section 4.2.4, however, this may not be the case for planets older than about 6 Gyr (Frank et al., 2014; Mojzsis, 2021). Any M dwarf still part of the distribution has an age much older than the Sun and even older than ∼6 Gyr, indicating that any geological activity on planets orbiting these old stars may already have vanished millions or even billions of years ago, at least for cosmochemically Earth-like planets.

OPEN IN VIEWER

FIG. 14.

Upper panel: The number of remaining GHZ stars that are above the lower limit, α at XUV , for our nominal case with as a function of birth age and separated into spectral classes. Here, cutoff birth ages can be seen for each spectral class (i.e., 0.05, 1.1, 2.7, and 11.5 Gyr for F, G, K, and M stars, respectively^a). The green dashed line shows the approximate age when geological activity at cosmochemically Earth-like planets will cease—no M dwarfs are younger than that age. Lower panel: The number of remaining GHZ stars for our nominal, minimum, and maximum cases as a function of stellar mass. Cutoff masses can be seen for each case. ^aNote that the temporal resolution of our model is 0.05 Gyr.

The lower panel of Figure 14 illustrates nominal, minimum, and maximum cases for x CO 2 , max = 10 % and shows the number of remaining stars as a function of stellar mass. As can be seen, the lower limit predominantly affects low-mass stars.

The numbers presented above, however, are for N₂-O₂-dominated atmospheres with x CO 2 , max = 10 % . If we set the mixing ratio to x CO 2 , max = 1 % instead, and implement the corresponding outer HZCL boundary from Schwieterman et al. (2019a) for pCO₂ = 0.01 bar, our numbers will decrease significantly. For this, we assume F X and F XUV threshold values in agreement with the upper atmosphere simulations by Johnstone et al. (2021a) for N₂-dominated atmospheres with x CO 2 , max = 1 % . As can be seen in Figure 35, N₂-atmospheres with 1% CO₂ will strongly expand already for F X ∼ 2.5 erg s − 1 cm − 2 . We may, therefore, define for our nominal and for our minimum cases. For the maximum case, we take to again account for any potential parameters that may enhance the stability of the atmosphere. However, since we have to assume a smaller HZCL with its outer boundary being defined through pCO₂ = 0.01 bar for all our cases, the number of remaining stars will be significantly smaller even for our maximum case. For the Sun at 5 Gyr, for instance, the mean HZCL distance of the maximum case will be located at only d 〈 HZCL 〉 = 1.05 AU, which clearly illustrates its small HZCL width.

By applying this criterion, only 0.73 ( + 2.54 / − 0.58 ) × 10 9 stars still remain that can host EHs covered by N₂-O₂-dominated atmospheres with x CO 2 , max = 1 % , a stellar fraction that accounts to 0.90(+3.13/−0.76)% of all stars in the galactic disk. The respective minimum stellar masses decrease toward M ⋆ = 0.65 ( + 0.09 / − 0.18 ) M ⊙ , which implies that no M dwarfs remain within the sample. Of all remaining stars, 2.90–58.32% are K dwarfs, while 33.49–79.17% and 8.19–17.92% are G and F stars, respectively. Figure 15 presents the nominal case for x CO 2 , max = 1 % .

OPEN IN VIEWER

FIG. 15.

Our nominal case for an N₂-O₂-dominated atmosphere with x CO 2 , max = 1 % and for F XUV , max = 15.0 erg s − 1 cm − 2 (upper and middle panels). The upper panel shows the distribution of remaining stars as a function of birth age and stellar mass, while the middle panel illustrates the distribution of stars within the different spectral classes as a function of birth age. The green dashed line again illustrates the approximate age at which geological activity will cease at cosmochemically Earth-like planets. The lower panel shows the number of remaining stars for all cases with x CO 2 , max = 1 % as a function of stellar mass. No M dwarfs are part of the remaining stellar distributions.

We emphasize that the actual number of stars suitable to host planets with N₂-O₂-dominated atmospheres will likely be lower, independent of the maximum level of pCO₂ chosen. Our calculation of the lower limit neither includes the effect of superflares nor nonthermal escape, the latter of which is triggered through a star’s stellar wind and CMEs. Since these processes are important particularly around M to late K dwarfs, these will further reduce the number of small-mass stars being part of the remaining stellar distribution. EHs around M- and late K-type dwarfs should, therefore, be exceedingly rare.