Potential long-term habitable conditions on planets with primordial H–He atmospheres

Abstract

Cold super-Earths that retain their primordial, H–He-dominated atmosphere could have surfaces that are warm enough to host liquid water. This would be due to the collision-induced absorption of infrared light by hydrogen, which increases with pressure. However, the long-term potential for habitability of such planets has not been explored yet. Here we investigate the duration of this potential exotic habitability by simulating planets of different core masses, envelope masses and semi-major axes. We find that terrestrial and super-Earth planets with masses of ~1–10 M_⊕ can maintain temperate surface conditions up to 5–8 Gyr at radial distances larger than ~2 au. The required envelope masses are ~10⁻⁴ M_⊕ (which is 2 orders of magnitude more massive than Earth’s) but can be an order of magnitude smaller (when close-in) or larger (when far out). This result suggests that the concept of planetary habitability should be revisited and made more inclusive with respect to the classical definition.

Main

The emergence of life as we know it is generally understood to require three building blocks: an energy source, access to nutrition and the presence of liquid water^1,2. Apart from these requirements, it remains unknown to what extent conditions on other planets need to resemble Earth to be habitable. Since the possibility of in situ solar system exploration, habitable conditions are considered on planets and moons that have liquid water oceans underneath a layer of ice³. These habitats would have to be inhabited by organisms that are adapted to high pressures and use a metabolism independent of photosynthesis. Such organisms have already been found on Earth; while today the majority of Earth’s biomass is concentrated on its surface (due to complex photosynthetic organisms such as land plants⁴), the subsurface biomass probably outweighed the surface biomass for most of geologic history⁵. Life has adapted to many other relatively extreme environments, such as the depths of the ocean at kbar pressures⁶, although it is unknown how many of these organisms live independently from life on the surface. In the search for life on extraterrestrial planets, it needs to be considered that life might manifest and thrive under conditions that would be considered extreme on Earth.

Observations of exoplanets have shown that the Solar System might not be a standard planetary system. This was initially implied after the first detections of hot Jupiters⁷ (although these turned out later to be rather rare but seemingly prevalent due to a detection bias) and further confirmed by the fact that super-Earths, absent in the Solar System, are very common^8,9,10. To date, super-Earths have only been observed at short orbits, but planet-formation models indicate that they are also prevalent at large radial distances¹¹, where they could retain a non-negligible H–He envelope, which is the original atmosphere accreted from the protoplanetary disk^12,13. At several au away from the host star, stellar radiation is insufficient to thermally evaporate a substantial amount of the atmospheric H–He¹⁴.

The primordial atmosphere, dominated by hydrogen and helium, would have insufficient greenhouse gases that are important on Earth, such as CO₂ or methane. However, if the atmosphere is massive enough, H₂ will act as a greenhouse gas. At sufficient pressures, the H₂ molecules undergo enough collisions to create a dipole moment, causing them to absorb the infrared radiation coming from the planet; this is known as ‘collision-induced absorption’^15,16. It could raise the surface temperature enough to allow for a liquid water ocean^17,18,19.

Past analytical work has shown that in the case of a massive atmosphere an intrinsic heat source is sufficient to warm the planetary surface, possibly even enabling unbound planets to be temperate¹⁷. Alternatively, planets with a smaller envelope can have stellar irradiation as the sole energy source, assuming that the accreted envelope is neither too massive nor too small for the given equilibrium temperature¹⁸. Recently, Madhusudhan et al.²⁰ investigated the mass–radius relationship of ‘Hycean worlds’, planets with a liquid water layer above a rocky core and below a H–He primordial atmosphere. Similar to previous studies^18,21, they found that liquid water can persist at a wide range of equilibrium temperatures. Therefore, the range of orbital distances where a liquid surface water can exist (the so-called habitable zone^22,23,24) is wider than in the case of Earth-like planets.

These past studies have concentrated on the static properties necessary for liquid water. However, liquid water conditions on these cold planets could be a transient phase¹⁹. It is important to consider atmospheric escape, as this can give limits to the separations at which H–He envelopes allow the existence and stability of liquid water. Here we study the temporal evolution in particular, to address the duration of a primordial H–He habitat.

We investigate how long surface liquid water could be permitted using an evolution model and define the duration of adequate surface pressure–temperature conditions as τ_lqw. First, we assume that the core mass, envelope mass and semi-major axis remain constant in time. We run a large grid of evolution models and vary these three parameters. During the planetary evolution, the equilibrium temperature changes due to the evolution of the host star, which is modelled as Sun-like²⁵. The planet’s intrinsic luminosity declines as the envelope and interior cool down and radiogenic components decay. For a given core mass we run 1,144 models, with the semi-major axis varied from 1 au to 100 au and the envelope mass from 10⁻⁶ to 10^−1.8 M_⊕.

Figure 1 shows these grids for models with core masses of 1.5, 3 and 8 M_⊕. Every grid point represents a planet with a given core mass, envelope mass and semi-major axis. The evolution models are run for 8 Gyr (see Methods for discussion). Coloured, circular grid points represent τ_lqw of at least 10 Myr. Planets with the longest duration, τ_lqw > 5 Gyr, are indicated in yellow. Within the range of considered core masses, the values of τ_lqw have the same trends with respect to envelope masses and semi-major axes. At orbits within ~10 au, long-term τ_lqw, indicated by the yellow data points, are distributed over envelope masses of 10⁻⁴ to 10⁻⁶. For envelopes smaller than ~10⁻⁵ the surface temperatures are temperate, similar to Earth’s. The stellar radiation is expected to be the dominant source of energy and there is a scaling relation between the received flux (through the semi-major axis) and the envelope mass for long-term τ_lqw. Beyond ~10 au, liquid water conditions follow the cooling of the interior. Planets with small envelope masses have liquid water conditions relatively early, while planets with more massive envelopes reach liquid water conditions later in their evolution.

Fig. 1: Duration of liquid water conditions for planets at a wide range of semi-major axes (1 au to 100 au) and envelope masses (10^−1.8 to 10⁻⁶ M_⊕).

Planets receive insolation based on the luminosity evolution of a Sun-like star. a–c, Core masses of 1.5 (a), 3 (b) and 8 M_⊕ (c). The duration of the total evolution is 8 Gyr. The colour of a grid point indicates how long there were continuous surface pressures and temperatures allowing liquid water, τ_lqw. These range from 10 Myr (purple) to over 5 Gyr (yellow). Grey crosses correspond to cases with no liquid water conditions lasting longer than 10 Myr. Atmospheric loss is not considered in these simulations. d, Results for planets with a core mass of 3 M_⊕, but including the constraint that the surface temperature must remain between 270 and 400 K. Every panel contains an ‘unbound’ case where the distance is set to 10⁶ au and solar insolation has become negligible.

Source data

In addition to τ_lqw, we show a case where an additional constraint is imposed: the surface temperature must remain below 400 K. This is based on the temperature limit for life on Earth²⁶. Figure 1d shows the results for this case for a planet with a 3 M_⊕ core.

Including atmospheric loss

The results presented so far do not include the effect of atmospheric loss. This is important for small planets close to their host stars, since they are much more vulnerable to stellar radiation^14,18,19. We apply the evaporation model discussed in the Methods to study which planets in the previous result would be stable against thermal atmosphere loss.

Figure 2 shows how the envelope mass changes with time with hydrodynamic escape (top) and Jeans loss (bottom) models. The Jeans escape model only has a notable effect on the mass of envelopes around small core masses of 1.5 M_⊕. For the Jeans model, we show the highest considered exosphere temperature of 2,000 K. This is relatively high, as exosphere temperatures for H–He compositions are estimated at ~1,000 K (ref. ²⁷). Since a higher exosphere temperature leads to a higher mass-loss rate, we consider this as an upper limit to Jeans escape. Lower values of 300 K and 1,000 K were also considered, but we find almost no effect on the envelope’s mass evolution. Even for the cases in which Jeans escape removes some of the planetary envelope we still find a similar τ_lqw, suggesting that the results presented in Fig. 1 are robust if escape occurs in the Jeans regime. The extreme ultraviolet (XUV)-driven hydrodynamic model, however, can substantially affect the evolution of small close-in planets over a wide range of parameters, as shown in Fig. 2.

Fig. 2: The evolution of the mass of the simulated envelope for two different atmosphere loss models: hydrodynamic escape and Jeans escape.

The top row shows the results of the hydrodynamic escape model and the bottom row of the Jeans escape model. The default values of the model are a core mass of 3 M_⊕ and initial envelope mass of 10⁻⁴ M_⊕ (M_{env, 0}) at a distance of 6 au. The left column shows how the envelope loss is affected by different semi-major axes, the middle column shows the effect of different envelope masses and the right column shows the effect of different core masses. Dashed lines show when a layer of water between the envelope and core would not be in the liquid phase; solid lines show when it would. Hydrodynamical escape can notably reduce or completely evaporate the envelope of some planets, which influences the duration of liquid water conditions. The Jeans escape model, even when a high exosphere temperature of 2,000 K is chosen, shows only a somewhat notable loss for a small core mass. of 1.5 M_⊕.

Source data

The results when the hydrodynamic escape model is included are shown in Fig. 3. In this case, we find that there are no long-term liquid water conditions possible on planets with a primordial atmosphere within 2 au. Madhusudhan et al.²⁰ found that for planets around Sun-like stars, liquid water conditions are allowed at a distance of ~1 au. We find that the pressures required for liquid water conditions between 1 and 2 au are too low to be resistant against atmospheric escape, assuming that the planet does not migrate at a late evolutionary stage.

Fig. 3: Duration of liquid water conditions for planets when including the effect of hydrodynamic atmospheric photo-evaporation.

The y axis shows the envelope mass after 8 Gyr. Escape inhibits liquid water conditions by removing the atmosphere for close-in planets with low initial envelope masses. Lower core masses are more affected.

Source data

Unbound planetary mass objects

The cases labelled as ‘unbound’ shown in Figs. 1 and 3 correspond to planets that are detached from their host star. Planets can be ejected after formation by gravitational interaction²⁸. Such unbound planets might be very common²⁹. Without the influence of a host star, these planets are not affected by post main-sequence evolution, which might terminate the duration of habitable conditions for bound planets. Instead, the habitable conditions would end when an internal heat source can no longer provide enough energy. Note that unbound planets are no longer protected by the stellar heliosphere and therefore galactic cosmic rays and γ rays could affect the planet’s structure and its potential habitability³⁰.

We simulate such cold planets on a longer timescale and investigate whether τ_lqw persists beyond the lifetime of a Sun-like star. The results are shown in Fig. 4. The simulations are performed up to 100 Gyr, when most of the radioactive nuclei have decayed and the intrinsic luminosity is nearly zero. The different results are caused by the longer integration time than for the unbound cases in Figs. 1 and 3.

Fig. 4: Duration of liquid water conditions on unbound planets.

The planets are simulated with different core masses and envelope masses. Grid point colours indicate how long there were continuous surface pressures and temperatures allowing liquid water, τlqw. The longest duration simulated was 84 billion years for a 10 M_⊕ core and a 0.01 M_⊕ envelope. If planets with a primordial atmosphere can host liquid water, the duration can be much longer on unbound planets since the internal heat source can evolve slower than the host star. Contour lines indicate the start of liquid water conditions for planets with τ_lqw > 100 Myr.

Source data

We find that many unbound planets are too hot shortly after their formation to host liquid water but do have the right conditions at later times and for long τ_lqw. In this case, the core mass is an important factor for τ_lqw, since we assume that the radiogenic heat source scales with core mass. The envelope’s mass influences the planetary cooling efficiency and therefore has a smaller, but non-negligible, effect on τ_lqw. Planets with cores more massive than 5 M_⊕ can have liquid water conditions lasting over 50 Gyr for envelopes with masses of ~0.01 M_⊕. The longest duration for liquid water is found to be 84 Gyr for a planet with 10 M_⊕ core and 0.01 M_⊕ envelope. The current set-up of our intrinsic luminosity model limits τ_lqw of unbound planets by the half-life of the radioactive components in the core, assuming a chondritic abundance. A potential planet with radioactive components with a longer half-life could therefore have an even longer τ_lqw. Many of the inferred τ_lqw values for unbound planets are much longer than the age of the universe. It should be noted, however, that it takes ~10 Gyr of cooling before these conditions are reached. Consequently, it could be that these types of planet are currently too hot for liquid water but will not be in a far future.

Our results indicate that habitable conditions underneath H–He-dominated atmospheres might be long-lasting and prevalent. Planets that retain their primordial atmosphere should have a relatively simple formation and evolution path compared with those with a secondary atmosphere. However, this is based on a lot of assumptions and simplifications on their formation and evolution. It is commonly assumed that habitable planets require a negative climate feedback to remain persistently habitable^22,23,31. However, this work suggests that the right set of initial conditions for H₂-rich super-Earths can alone lead to persistently temperate surface conditions over Gyr timescales. There are no accurate predictions on the occurrence of super-Earth-sized planets with these initial conditions, but it is likely enough that these alternatively habitable planets constitute a fraction of the habitable worlds in the galaxy.

We made the assumption that the combination of core masses and envelope masses of these initial conditions can form. Theory³² and simulations^11,33 indicate more massive envelopes are likely to form beyond the snow line. Studies of observed close-in super-Earths suggest that these planets consist of H–He envelopes of several percent in mass^34,35. Rogers and Owen¹³ showed that the inferred post-formation H–He mass fraction peaks at ~4%, but extends well into the range needed for long τ_lqw (envelope mass M_env ≤ 10^−2.5 M_⊕). The same work showed that formation models predict too-high envelope mass fractions, motivating the improvement of planet-formation models. Less-massive envelopes could be the result of additional mechanisms, such as core-powered mass loss^36,37 or collisions^38,39. We suggest that future studies should investigate whether the planetary evolution of unbound planets differs from that of bound planets and if so how it may affect their potential habitability. A detailed study on the formation likelihood of planets with liquid water conditions underneath H–He envelopes that are observable at the current time is beyond the scope of this research and we hope to address it in future studies.

Future work should investigate in detail the formation likelihood of planets with the right initial conditions. The possibility of their formation around M-dwarfs should also be considered since these stars have a prolonged pre-main sequence with higher UV fluxes. Planets around such stars are expected to lose more of their primordial envelope and water^40,41. Therefore, future work should consider the initial volatile abundances of such planets to predict the water mass fraction at different evolutionary stages^42,43. In addition, while a sufficient amount of water left after atmospheric escape is required, the amount of water cannot be too high, to prevent the formation of ice layers between the silicate core and the liquid water. The exchange of nutrients would then be inhibited^44,45 or at least limited in the possible scenario that the ice layer is fully convective⁴⁶.

Another simplification is that the planets accrete a solar-composition envelope and retain this composition. However, processes during evolution can increase metallicity. This can lead to a different total greenhouse effect, with contributions of CO₂ as well as collision-induced absorption by H₂–N₂ and H₂–CO₂ (refs. ^40,47,48,49). We have explored the sensitivity of our result to different compositions in the Supplementary Information.

Life on the type of planet described in this work would live under considerably different conditions than most life on Earth. The surface pressures in our results are on the order of 100–1,000 bar, the pressure range of oceanic floors and trenches. There is no theoretical pressure limit on life, and some of the most extreme examples in Earth’s biosphere thrive at ~500 bar (ref. ⁵⁰). These habitats also receive a negligible amount of direct sunlight, and therefore photosynthesis would not be an optional mechanism to provide for metabolism. Chemoautotrophic life on Earth⁵¹ would be a more likely analogue to possible life on this type of planet. Earth is inhabited by lifeforms that might be adapted to life underneath a H–He envelope⁵². Since chemotrophic life arose before photoautotrophic lifeforms, it can exist independently as it presumably did on Earth⁵³. However, it is not trivial if the emergence of life could happen on such planets in a similar way as on Earth, given that the planets we simulated spend a longer time of their evolution being too hot⁵⁴. Furthermore, the advent of photosynthesis on Earth introduced a much more productive metabolism not dependent on pre-existing chemical energy gradients. Consequently, planets dominated by photosynthetic life probably produce more readily observed signatures of their inhabitation. The lack of sunlight means there would also be no temporal variation from a day–night cycle or season. This would lead to stable conditions (which could influence the kind of lifeforms that are best adapted^2,55). Another property that can affect planetary habitability is the presence of a magnetic field⁵⁶. It is yet to be determined whether the planets we consider in this study are expected to generate a dynamo, and we hope to address this in future research.

More recent interpretations of planetary habitability do not only correspond to harboring life, but also to the ability to detect it from Earth⁵⁷. Currently it is still challenging to observe small planets at large radial distances, but great advancement is expected from future telescope missions. The relatively large-scale height of H–He atmospheres could make it possible for instruments such as the James Webb Space Telescope and Ariel or the Extremely Large Telescope to detect and characterize biomarkers in such atmospheres⁵². While life in a H₂-dominated atmosphere environment could produce biomarkers⁵⁸, one challenge is that habitats that lack photosynthetic life do not produce a chemical disequilibrium, but rather destroy it by their metabolism⁵⁹. Future studies should predict if certain chemical disequilibria are biomarkers or anti-biomarkers for these specific habitats. Because of the importance of the system’s temporal evolution, age determination by PLAnetary Transits and Oscillations of stars (PLATO)⁶⁰ is also critical for the assessment of planetary habitability. The Roman space telescope, using gravitational microlensing, could detect colder exoplanets and even unbound planets^61,62 to constrain their population. We therefore expect that our understanding of this exoplanetary population and its potential habitability will substantially improve in the near future.

Methods

We simulate the gaseous envelope using a spherically symmetric model that solves the structure equations of a planet’s interior assuming hydrostatic equilibrium^63,64,65:

$$\frac{{{{\mathrm{d}}}}m}{{{{\mathrm{d}}}}r}=4\uppi {r}^{2}\rho$$

(1)

$$\frac{{{{\rm{d}}}}P}{{{{\rm{d}}}}r}=-\frac{Gm}{{r}^{2}}\rho$$

(2)

$$\frac{{{{\rm{d}}}}\tau }{{{{\rm{d}}}}r}={\kappa }_{{{{\rm{th}}}}}\rho$$

(3)

$$\frac{{{{\rm{d}}}}l(r)}{{{{\rm{d}}}}r}=0.$$

(4)

This relates the local cumulative mass m to the radius r within the planet. P and ρ are the local pressure and density, respectively, and G is the gravitational constant. κ_th is the Rossland mean opacity. l is the local intrinsic luminosity: the energy flux going through a shell of radius r. Contributions to the luminosity due to contraction and cooling of both the envelope and the iron–silicate core (Supplementary Fig. 1) are included, as well as radiogenic heating (Linder et al.⁶⁵). The implication of equation (4) is that the luminosity does not vary through the gaseous envelope. This simplification is valid when the envelope has a small contribution to the total luminosity compared with the core (Intrinsic luminosity).

The temperature (T) gradient at high optical depth, τ, is represented either by the adiabatic or radiative gradient, using the Schwarzschild criterion^66,67. In the case of a convective region, we use the adiabatic gradient:

$$\frac{{\mathrm{d}}\,T}{{\mathrm{d}}\,P}=\frac{T}{P}{\left(\frac{{{\mbox{ln}}}T}{{{\mbox{ln}}}P}\right)}_{{{{\rm{adi}}}}}.$$

(5)

At low optical depths, a double-grey atmosphere model⁶⁸ is adapted. This atmosphere model uses the assumption that the spectral energy distribution of received stellar radiation is mostly optical light, while the outgoing radiation of the planet is mostly in the infrared. The temperature as a function of optical depth is given by:

$$\begin{array}{l}T{(\tau )}^{4}=\frac{3{T}_{{{{\mathrm{int}}}}}^{4}}{4}\left\{\frac{2}{3}+\tau \right\}+\frac{3{T}_{{{{\mathrm{eq}}}}}^{4}}{4}\left\{\frac{2}{3}+\frac{2}{3\gamma }\left[1+\left(\frac{\gamma \tau }{2}-1\right){{\mathrm{e}}}^{-\gamma \tau }\right]\right.\\\left.\qquad\quad+\frac{2\gamma }{3}\left(1-\frac{{\tau }^{2}}{2}\right){E}_{2}(\gamma \tau )\right\}\end{array}$$

(6)

The equilibrium temperature T_eq is defined as ${T}_{{{{\mathrm{eq}}}}}={T}_{\star }\sqrt{\frac{{R}_{\star }}{2a}}{\left(1-{A}_{{{{\mathrm{B}}}}}\right)}^{\frac{1}{4}}$, where ${T}_{\star }$ and ${R}_{\star }$ are the effective temperature and radius of the star, respectively, a is the distance between planet and star and A_B is the Bond albedo fixed to Jupiter’s value of 0.343. The intrinsic heat, T_int, is defined as ${T}_{{{{\mathrm{int}}}}}=\left(\frac{{L}_{{{{\mathrm{int}}}}}}{{\sigma }_{{\mathrm{B}}}4\uppi {r}^{2}}\right)$. Here L_int is the total intrinsic luminosity and σ_B is the Stefan-Boltzmann constant. E₂(γτ) is a form of the exponential integral ${E}_{n}(z)\equiv \int\nolimits_{1}^{\infty }{t}^{-n}{{\mathrm{e}}}^{-zt}{\mathrm{d}}t$, where n = 2, z = γτ and t is the integration variable. The greenhouse effect is parametrized as the ratio between the optical opacity and the infrared opacity as γ = κ_vis/κ_th. The level of greenhouse effect is thus contained in the parameter γ. We use a parameterization of γ (table 2 in Jin et al.⁶⁴). Since most of our planets have equilibrium temperatures that are lower than the low limit in this table (T_eq,min = 260 K), we extrapolate the table for lower γ values. It should be noted that the term in equation (6) containing the intrinsic heat becomes the dominant term in the extrapolated regime of low equilibrium temperatures, and therefore different γ values have a small effect on the resulting temperature. We therefore find that the exact value of γ in the extrapolated regime has a negligible effect on the inferred state of water underneath the atmosphere. The atmosphere then reduces to a classical Eddington grey atmosphere. At high optical depths, the temperature in convectively stable regions is calculated with diffusive radiative energy transport as:

$$\frac{{{{\mathrm{d}}}}T}{{{{\mathrm{d}}}}r}=-\frac{3{\kappa }_{{\mathrm{R}}}\rho l}{64\uppi {\sigma }_{{{{\mathrm{B}}}}}{T}^{3}{r}^{2}}.$$

(7)

The transition between the low-optical-depth regime (equation (6)) and high-optical-depth regime (equation (7)) is determined by two boundaries, between which we interpolate the temperature. These two boundaries are given by: $\tau =10/\sqrt{3}\gamma$ and $\tau =100/\sqrt{3}\gamma$. The transitions occur at 10 bar and 100 bar, which means that planets that have an envelope of ≾10⁻⁵ M_⊕ are completely in the low-optical-depth regime.

The model uses Rossland-mean opacities for a solar- or a scaled-solar-composition gas. One source of opacities is the molecular, grain-free tables⁶⁹. These are tables in temperature–density space and depend on the metallicity [M/H], covering a range of [M/H] = 0–1.7 (close to 50 times solar). They include collision-induced absorption, which is expected to dominate at low temperatures. The lower limit on temperature for this opacity table is 75 K. Our simulations sometimes contain planets with lower temperatures in certain regions of the atmosphere (down to 50 K). In this case, we extrapolate the opacities as ∝T², as suggested by the temperature dependency close to 75 K. To assess the consequences of the extrapolation, we also repeated our simulations with a fixed opacity in the extrapolated regime. This made no difference to our final result, since the planets that we find to allow liquid water do not have a substantial part of their atmosphere in the extrapolated regime.

We use a non-ideal equation of state⁷⁰ for H–He¹⁵. For H₂O we use AQUA⁷¹.

Solid core

Our internal structure model⁷² solves equations (1) and (2) for the solid core. We assume for all cores a silicate to iron ratio of 2:1, similar to Earth, and assume that there are no ices. It seems likely that planets forming beyond the snow line accrete a substantial amount of water ice. On the other hand, several mechanisms dehydrate the planetary building blocks⁷³. Including ices in our model has several possible consequences. First, it changes the mass–radius relationship. Since the radii of water-rich planets are more sensitive to temperature, we expect that the core radii will vary more during the evolution. Second, at a fixed total core mass, the presence of water would reduce the silicate mantle mass and thus the amount of estimated radiogenic heating (Variations in intrinsic luminosity). Lastly, it would allow more energy to be stored in the core, which is released over longer timescales due to the fact that ice has a higher specific heat capacity than silicates. We plan to account for different planetary compositions and investigate their effect on planetary evolution in future research.

The core model uses a modified polytropic equation of state⁷⁴ that relates the density ρ and the pressure P as:

$$\rho (P)={\rho }_{0}+c{P}^{n}$$

(8)

ρ₀, c and n are parameters specific to the material⁷⁴. We use MgSiO₃ for silicates⁷². The external pressure that the envelope imposes on the core is included, although this does not have a notable effect on our calculations as the most massive envelopes considered are only 0.01 M_⊕.

Intrinsic luminosity

The intrinsic luminosity of the planet and its temporal evolution are included in the simulation. We estimate the initial value of this intrinsic luminosity and its evolution in the following way.

First, we estimate the initial intrinsic luminosity that the planet has shortly after formation. We use an analytical fit⁷⁵, which uses the results of numerous planet-formation simulations. The intrinsic luminosity of planets in the formation simulation is calculated by solving the one-dimensional structure equations, including the heating from cooling and contractions as well as accretion of gas and planetesimals. The analytical approximation yields the intrinsic luminosity as it depends on core mass, envelope mass and age of the planet. We fix the starting age of the planet at the same value as the start of our simulation, namely 20 Myr. The initial luminosity is then found by interpolating in core and envelope mass. Our simulations also show how much of the luminosity is generated in the core or the envelope (an example of which is shown in Fig. 4). This has also been described in detail in Linder et al.⁶⁵. It is found that the luminosity contribution of the H–He envelope is much smaller than the contribution of the solid core, in agreement with earlier work^34,65. This justifies the assumption of a uniform luminosity in the envelope made by equation (4). If a substantial part of the total luminosity comes from envelope cooling and contraction, this assumption would not be valid.

The second contribution used to calculate the total intrinsic luminosity is based on radiogenic heating⁷². We model the heat due to the decay of radioactive nuclides by assuming an initial abundance of radioactive material (⁴⁰K, ²³²Th and ²³⁸U) of a chondritic composition. Then we assume that the total abundance of the radiogenic heating scales with the mass of the silicate mantle. The radiogenic luminosity becomes:

$${L}_{{{{\rm{radio}}}}}(t)={L}_{{{{\mathrm{radio}}}},\oplus }(t)\times \frac{{M}_{{{{\mathrm{mantle}}}}}}{{M}_{{{{\mathrm{mantle}}}},\oplus }}.$$

(9)

Here M_mantle,⨁ and L_radio,⨁(t) are Earth’s mantle mass and Earth’s radiogenic luminosity at time t, respectively. M_mantle is the mantle mass of the modelled planet, which we assume to be 2/3 of the core mass. Once the initial luminosity is found, its temporal evolution is given by energy conservation, that is, by the condition that the total energy difference between two points in time is equivalent to the luminosity radiated between them⁶³. The radioactive component evolves according to the decay time of the radioactive nuclides.

Figure 1 shows an example of the evolution of the intrinsic luminosity of a 3 M_⊕ core and a 0.001 M_⊕ envelope for a planet that receives a negligible amount of stellar radiation. After 300 Myr, the radiogenic luminosity becomes the dominant source of intrinsic heat. Since the abundance of radioactive nuclides in other planets is unknown⁷⁶, we also consider cases where these are a factor of ten larger/smaller, which are presented in the Supplementary Information under Variations in intrinsic luminosity.

Supplementary Fig. 2 shows how the atmosphere structure evolves in 3 planets from 50 Myr to 5 Gyr. All planets have a core mass of 3 M_⊕ and an envelope of 10⁻³ M_⊕. No atmosphere-escape model is implemented; the evolution is due to the host star and the intrinsic luminosity. The planets at 1 and 10 au remain too hot for liquid water conditions, while the planet at 100 au reached liquid water conditions after 1 Gyr.

Atmosphere loss

We model atmospheric loss through two distinct thermal escape models: one based on Jeans escape and one based on hydrodynamic escape.

Jeans escape is applicable for planets that remain in hydrostatic equilibrium everywhere^27,77,78,79. The fraction of escaping particles depends on the escape velocity, particle number density and the temperature (T_exo) at the exobase. For the escape velocity, we use the total radius and total mass of the planet and therefore assume the escape velocity at the exobase is the same as at the outer radius. We also assume a 100% hydrogen composition for the number density and neglect that helium particles are heavier and thus have lower velocities at the same temperature when assuming hydrostatic equilibrium. This results in an overestimate of the escape rate, but using an upper limit of Jeans escape serves our goal of investigating whether Jeans escape has an important effect on liquid water duration.

The exosphere temperature is complex to estimate, as it is determined by the absorption of X-rays and XUV and by the composition at the exosphere⁸⁰. H–He-dominated atmospheres are expected to have a relatively warm exosphere of ~1,000 K (ref. ²⁷), while planets dominated by, for example, CO₂, H₂O or N₂ should have a colder T_exo of ~300 K (ref. ⁸¹). We remain agnostic about the specific exosphere temperature. Instead, we put the exosphere temperature to a fixed value during the simulation (T_exo,min) and use a range (300–2,000 K). If the temperature of the outer radius of the atmosphere (where τ = 2/3) is warmer than the assumed exosphere temperature, we allow the exosphere temperature to be increased to this value at any moment in the evolution. However, we find this never to be the case.

The second escape model is based on hydrodynamic escape^64,82. Hydrodynamic escape is mostly driven by X-ray radiation at early times and later by XUV radiation⁸³. We assume X-ray-dominated evaporation until a threshold of XUV flux is reached^64,83. The evolution of the stellar X-ray and XUV luminosity are from solar models⁸⁴. During X-ray-dominated evaporation, the loss is energy limited⁸⁵. We estimate the mass loss by^64,86:

$${\dot{M}}_{{{{\mathrm{X-ray}}}}}=\epsilon \frac{\uppi {F}_{{{{\mathrm{X-ray}}}}}{R}_{2/3}^{3}}{G{M}_{{{{\mathrm{pl}}}}}\,{K}_{{{{\mathrm{tide}}}}}}$$

(10)

where F_X-ray is the received X-ray flux, ${R}_{2/3}^{3}$ is the radius at optical depth of 2/3, M_pl is the planet’s total mass and K_tide is an extra factor to account for a higher mass loss when planets have their Roche-lobe boundary close to the surface⁸⁷. Since we simulate relatively small planets at large distances, K_tide has a negligible effect. In contrast to earlier work^64,82, we now parameterize the efficiency factor ϵ by the escape velocity⁹.

The XUV-dominated loss can also be energy limited. In this case the mass loss is also calculated with equation (10), but using the XUV flux and the radius where the optical depth in UV equals 1 (ref. ⁸⁸). When the XUV radiation is high, a substantial amount of the heat can be lost to radiation. In this case, the radiation/recombination limited (rr-limited) description⁸⁸ is used. Our model computes both the energy-limited and rr-limited mass loss and applies the lowest value.

Liquid water conditions

During the simulations, we determine whether there are liquid water conditions by checking the pressure and temperature at the bottom of our simulated atmosphere. We compare these with the phase diagram of water to determine if a water layer could be in the liquid phase⁸⁹. We define the duration of liquid water conditions as the time that water is permitted in the liquid phase without interruptions, referred to as τ_lqw. In some of our results, we apply an extra constraint to the liquid water conditions: that the surface temperature remains below 400 K. This is based on the observation on Earth that terrestrial life thrives best at temperatures of around 300 K and an upper limit for the chemistry of life is estimated at ~400 K (ref. ²⁶).

Comparison with Pierrehumbert and Gaidos¹⁸

To test our model in time-independent calculations, we first compare our models with those of Pierrehumbert and Gaidos¹⁸ in Supplementary Fig. 3. We calculate the necessary atmosphere mass for a surface temperature of 280 K. In addition to our own model (herein referred to as ‘COMPLETO’) described in the sections above, we use another atmosphere model: PETITcode⁹⁰. PETITcode is also a one-dimensional radiative–convective equilibrium code, but with a wavelength-dependent treatment of radiative transport. It uses specific chemical abundances and assumes a chemical equilibrium. At low temperatures, however, non-equilibrium effects can play an important role in the atmosphere. It assumes a uniform gravity throughout the atmosphere.

The simulations of Pierrehumbert and Gaidos used a constant gravity of 1,700 cm s⁻², without an interior heat source. To match this, we set the core mass at 3 M_⊕, which resulted in a radius of 1.32 R_⊕. The intrinsic heat value was set to a negligible temperature (1 K internal temperature for PETITcode, 10⁻⁶ L_Jup, where L_Jup is Jupiter’s current luminosity at L_Jup = 3.35 × 10¹⁹ erg sec^–1, for COMPLETO). Furthermore, Pierrehumbert and Gaidos¹⁸ used a pure H₂ composition. In PETITcode we set the metallicity to 0.01 times solar, so that there is predominantly hydrogen and a non-negligible fraction of helium, though we do not expect this helium to influence the opacities. In COMPLETO we use solar composition since ref. ⁶⁹ only provides tabulated opacities for [M/H] = 0–1.7. The age of the host star is set at 5 Gyr.

The results of the PETITcode models agree very well with those of Pierrehumbert and Gaidos¹⁸. The results are also shown with the inclusion of an intrinsic heat source of 35 K, which is the average intrinsic temperature we find for a 3 M_⊕ planet when applying our intrinsic luminosity model. These results show that, as expected, at far semi-major axes this intrinsic heat source can reduce the amount of atmosphere needed to warm the surface. Quite good agreement is found inside of about 7 au. Outside, higher pressures are found, probably a consequence of extra absorption of incoming visual light in the solar-composition atmosphere relative to more H–He.

Illustrative cases of evolving models

Supplementary Fig. 4 shows the evolution of surface pressures and temperatures for different models, and if the pressure and temperature allow liquid water (solid lines) or not (dashed lines). We consider typical low-mass planets (M_core of 1.5 to 8 M_⊕) with envelopes of 10⁻⁵ to 10⁻³ M_⊕. The duration of continuous liquid water conditions, τ_lqw, is therefore the integration of the (longest) solid line.

Supplementary Fig. 4a shows the effect of changing the semi-major axis. Planets relatively close to their star at 2 au are too hot for the first 500 Myr, when the combination of high irradiation (from a small semi-major axis) and a high intrinsic luminosity (from a young age) lead to high temperatures. These planets develop liquid water conditions at a later stage, while the same planet further away (at 6 au or 10 au) can host an ocean earlier. The effect of changing the atmospheric mass is shown in Supplementary Fig. 4b, where the core mass and semi-major axis are fixed. This directly influences the pressure at the bottom of the atmosphere. An envelope of 10⁻³ M_⊕ leads to temperatures that are too high. A smaller envelope of 10⁻⁴ M_⊕ has liquid water conditions after 70 Myr, at high pressures and temperatures. An even smaller atmosphere of 10⁻⁵ M_⊕ has liquid water conditions at the beginning of the evolution, but then becomes too cold. In Supplementary Fig. 4c we show the effect of the core mass. Changing the core mass has two different effects on our model. First, it will result in different gravities in the atmosphere model with consequences on the radiative–convective profile. Second, the intrinsic luminosity depends on the core mass.

Dependence on model parameters

In this section, we test the sensitivity of our results to the assumptions we made. We show in Supplementary Fig. 5 how different model parameters influence the duration of liquid water conditions in comparison with the nominal case presented in Fig. 1.

The effect of including a 0.5 ice mass fraction in the core is shown in Supplementary Fig. 5b. The nominal case assumes no ice in the core. The biggest consequence of assuming the core is half made up of ice is that it leads to an estimated mantle that is half as massive and therefore an amount of estimated radiogenics that is a factor of two lower.

Variations in intrinsic luminosity

The second parameter we investigate is the intrinsic heat. The value and evolution of the intrinsic heat is an important, if not the main, contributor to habitable surface temperatures for the planets considered in this work. We change our intrinsic heat model in two separate ways: by changing the initial luminosity due to cooling and contraction at the start of our simulation and by changing the magnitude of the radiogenic luminosity.

A change of the initial luminosity (excluding the radiogenic component) has little influence on the long-term surface conditions. Planets with more (less) initial luminosity will start with more (less) entropy. This results in more (less) efficient cooling and contractions. In our simulations, it takes around 100 Myr for planets with different initial luminosities to converge to the same luminosity. Since we find τ_lqw on the timescale of billions of years, it remains unaffected even when we change the initial luminosity by a factor of four.

The radiogenic luminosity component is independent to the rest of the intrinsic luminosity sources. As we show in Supplementary Fig. 1, this component becomes the dominant term on long timescales and therefore we expect it to be critical in the determination of τ_lqw. We perform our simulations but multiply or divide the radiogenic heat source by a factor of 10. This would correspond with a planet that has 10 times more or less radioactive material as a fraction of the mantle core, respectively.

In Supplementary Fig. 5c,d we show how τ_lqw is influenced by changes in L_radio of a factor of 10. Up to ~10 au, τ_lqw seems mostly unaffected. The exception are planets with an envelope of ~10⁻⁵ at ~7 au that are too hot for liquid water in the nominal case but cool enough in the low L_radio case. At distances on the order of ~10–100 au, τ_lqw depends on L_int. Since the colder, more massive planets in our result are generally too hot for liquid water, the case of 0.1 × L_radio increases the duration for planets with a larger envelope. We only took into consideration the effect of a different magnitude of L_radio, which would correspond to a different total amount of radioactive particles. We did not take into account that a different composition of radioactive material could additionally lead to a different half-life.

Variations in composition

It is likely that planets that retain their primordial atmosphere still have varying degrees of metallicity compared with solar metallicities. Also, the relative abundances of elements will vary. Given the uncertainty in ranges of possible exoplanet compositions, we investigate how they potentially affect our results by changing two general parameters: the greenhouse parametrization γ in equation (6) and the opacities.

The ratio between visible and thermal opacities $\gamma =\frac{{\kappa }_{{{{\rm{vis}}}}}}{{\kappa }_{{{{\rm{th}}}}}}$ incorporates the greenhouse effect in the atmosphere model. Ideally, the value of γ takes into account condensation of any molecule that absorbs in the infrared, with contributions of all possibly relevant species such as CH₄ and CO₂.

In the nominal case, γ is extrapolated from table 2 in Jin et al.⁶⁴ from the effective temperature. Therefore, the value of γ changes in time. As discussed above, the atmosphere model has a small dependency on γ when the intrinsic temperature is substantially larger than the equilibrium temperature. This is the case for young planets and/or those that are at large semi-major axes.

Supplementary Fig. 5e,f shows the effect of fixing the value of γ to 0.01 or 0.001. It is in line with the double-grey atmosphere model that the value of γ is an important parameter for determining τ_lqw at shorter distances (>3 au), but does not affect the planets at larger distances.

Finally, we investigate the effect of enhanced infrared opacities in Supplementary Fig. 5g. We increase the calculated opacities by a factor of 10. Enhanced infrared opacities could be the result of relatively more greenhouse gases present in the envelope, such as CH₄, NH₃ or H₂O. Our results show that enhanced infrared opacities lead to larger surface pressures and therefore liquid water conditions would be found on planets with smaller envelopes. In reality, different greenhouse gasses would absorb at different wavelengths. To account for the contributions of specific enhanced greenhouse gasses would require a more complicated model.

Comparison with planet population synthesis

Given the uncertainty in planet formation, it is not possible to give an accurate prediction of the occurrence rate of planets that satisfy the right initial conditions leading to liquid water conditions. Nevertheless, we use the New Generation Planet Population Synthesis (NGPPS)^11,91 to compare the predicted initial conditions with those that favour long-term liquid water conditions. The NGPPS consists of 1,000 systems all formed around a 1 M_★ star and integrated to 100 Myr. After the integration time, 34,635 embryos are formed.

Most planets accrete envelopes that are too big for liquid water conditions, making a potential surface temperature too hot. Nevertheless, planets with the required initial conditions do occur. Out of the 34,635 embryos formed, 6,851 have envelope masses of 10⁻⁶ < M_env/M_⊕ < 0.01 and 5,160 have 10⁻⁶ < M_env/M_⊕ < 0.001. Out of these, most have a small core mass (below 1 M_⊕), with only 187 planets having a core mass between 1 and 10 M_⊕. There are several reasons why the envelope mass can be overestimated in the NGPPS. First of all, it applies one-dimensional, hydrostatic models. More realistic models find reduced gas accretion because of gas exchange^92,93. The assumed grain opacities can also be an underestimation, leading to too efficient gas accretion. Other assumptions on formation and evolution could also lead to smaller primordial envelopes than estimated in the NGPPS, for example, core-powered mass loss^36,37 or collisions^38,39.

The fact that the NGPPS predicts (too) high envelope masses is in line with the finding that planet formation theory overestimates envelope masses in comparison with observations¹³. It is therefore desirable to identify the missing physical/chemical processes that lead to the higher envelope masses. We also suggest that more sophisticated planet formation and structure models are required to estimate the occurrence rate of planets with liquid water at the current time.