© The Authors 2024

1 Introduction

According to the Planetary Habitability Laboratory¹ (PHL), to date there are 59 known rocky exoplanets that are orbiting in the habitable zones (HZs) of their host stars. There are different approaches to investigating the habitability of exoplanets. Since habitability on the surface is strongly dependent on the presence of an atmosphere, detecting and measuring the composition of exoplanetary atmospheres – which is in the scope of the recently launched James Webb Space Telescope (Greene et al. 2019; Birkmann et al. 2022) – is one of the most important tools for studying the habitability of exoplanets today. However, in the absence of an atmosphere, surface life is unlikely to develop (Lammer et al. 2009). Atmospheric loss processes can result in planets that are completely deprived of their atmospheres. These processes could be particularly significant for exoplanets orbiting close to M dwarf stars. Nevertheless, planets without a significan atmospheric layer can still provide the conditions necessary for the development of life. Such planets could be habitable given the presence of an ocean under the surface, as in the case of Enceladus. Therefore, it is important to explore other possible ways of studying the habitability of exoplanets that do not rely on the presence of an atmosphere.

Relatively little attention has been paid to the possibility that life might flourish in subsurface liquid water oceans on exoplanets. In the presence of water beneath a layer of surface ice, the interior of a planet might be an ideal place for life to develop, because the surface ice layer could provide protection against energetic particles (Paranicas et al. 2009; Pavlov et al. 2019). Since water can exist beneath the surface even in the absence of an atmosphere, investigating the internal structures of exoplanets may open a new window for habitability studies that does not depend on atmospheric considerations. Most of the currently known rocky exoplanets orbit close to their stars, which likely prevents them from retaining their water. Habitable zone exoplanets are ideal targets since it is possible to investigate both surface and subsurface habitability. High H₂O content, for example, may imply either subsurface water reservoirs or, in the presence of an Earth-like atmosphere, global surface oceans. It is important to emphasize that investigating the possibility of subsurface water reservoirs does not exclude the presence of planetary atmospheres or liquid oceans on the surface, but extends possible habitable regions to subsurface regimes. Furthermore, investigating the interior structure of exoplanets may contribute to improving our understanding of the relationship between planetary formation and composition as well as to helping constrain models of planetary evolution (Valencia et al. 2006; Elkins-Tanton & Seager 2008; Grasset et al. 2009).

However, determining the interior structure of an exoplanet is a significant challenge, because the only available relevant information is either the radius or the mass of the planet or in some cases both of these parameters. Uncertainties in the planetary masses and radii can result in a wide range of possible values for density. Considering that there is a wide variety of interior structures that correspond to the same density, these uncertainties in the planetary parameters result in an extended range of plausible scenarios for the internal composition (e.g., Rogers & Seager 2010; Unterborn et al. 2015; Barr et al. 2018; Dobos et al. 2019). Precise measurements of planetary masses and radii could narrow this range and improve our understanding of planetary interiors.

The habitability of planet interiors depends on the availability of liquid water. Icy moons in our Solar System, such as Enceladus or Europa, are examples of rocky bodies with an extended ocean underneath a thick layer of ice (Khurana et al. 1998; Zimmer et al. 2000; Iess et al. 2014; van Hoolst et al. 2016). In the case of these moons, the oceans remain liquid due to the heat flux provided by radiogenic heating and tidal dissipation (Nimmo et al. 2007). The recently launched ESA spacecraft, the Jupiter Icy Moons Explorer (Juice) is the next step towards studying the interior of icy bodies (Grasset et al. 2013). Juice will investigate the icy moons of Jupiter in order to gain a deeper understanding of the interior structure, formation, and potential habitability of these bodies. Analogous to these icy moons in our Solar System, H₂O rich exoplanets with strong internal heating may maintain a layer of liquid water under the planetary surface. It is important to note that a substantial H₂O-reservoir may indicate a surface ocean or ice layer, depending on the surface temperature of the planet. Even in the case of HZ exoplanets, in order to gain more insight into the habitability of these bodies, it is crucial to gain knowledge about the water content of the planets.

For close-in exoplanets with non-zero eccentricities, tidal dissipation is caused by the rising and lowering of a tidal bulge on the planet. The rising and lowering of the bulge during the course of the planet’s orbit is resisted by the planet’s own internal friction. This results in heat (Cassen et al. 1979). The time variability of the tidal forces exerted by the host star on the planet can be a significant contributor to the total energy budget of the planet and may play an important role in habitability as well (Ferraz-Mello et al. 2008; Jackson et al. 2008a,b; Běhounková et al. 2011).

Another process that may contribute to internal heating is the decay of long-lived radioactive isotopes. Radiogenic heating is a significant contributor to the internal energy budget of the Earth (Schubert et al. 2001). Gando et al. (2011) found the contribution of radiogenic heating to the heat flux of the Earth to be ${20.0}_{-8.6}^{+8.8}$ TW, which converts to a surface heat flux of around 0.04 W/m². They concluded that the decay of radioactive elements is responsible for approximately half of the total internal heat flux of the Earth.

The amount of radiogenic heating is strongly dependent on the planetary composition, in terms of the mass fraction of rock and the abundances of radioactive isotopes, chiefly ²³⁸U, ²³⁵U, ²³²Th, and ⁴⁰K. The radiogenic abundances in exoplanets may be quite different from what we see in our own Solar System, and could depend on the metallicities of the host star. Unterborn et al. (2015) found that the Th-abundance of solar analogues are generally higher than the solar value, suggesting that exoplanets orbiting these stars can experience higher heat fluxes from the decay of radioactive elements than planets in our Solar System. Similarly, Botelho et al. (2018) found that the Th abundance in the galactic thin disk was on average super-solar even for earlier times during the evolution of our Galaxy, and that the Sun is deficient of Th when compared to solar twin analogues.

Here, we present our studies on the internal structures, water content, and heat budgets of a selected sample of well-characterized exoplanets. In Sect. 2, we describe our approach for modelling the interior structures of exoplanets and our methods for computing tidal and radiogenic heating. In Sect. 3, we present our results and assess the possibility of the existence of liquid water reservoirs in the investigated HZ exoplanets. We summarise our main conclusion in Sect. 4.

2 Methods

2.1 Physical and orbital parameters

We studied HZ rocky exoplanets to assess their habitability based on their structures and internal heat production. The TRAPPIST-1 planets were ideal candidates for interior modeling because their radii and masses are known to a high precision. We also included other HZ rocky bodies in our research since they are of key importance in habitability studies (see Table 1 for the full list). Rocky exoplanets that are orbiting far from their host stars outside of the HZ could be considered for investigation from the aspect of subsurface water reservoirs. However, only a handful of such planets are known, and all of them were discovered by gravitational microlensing (e.g. Gould et al. 2014; Jung et al. 2019; Zang et al. 2021). Due to the small chances of further characterisation, these planets are not considered in this study. Although the number of presumed rocky HZ exoplanets is larger than the number of exoplanets investigated in this research, we only modelled those bodies that had a probability of over 5% of having a rocky composition (see Sect. 3.4). To determine this, we used the Forecaster model, which uses an empirical relationship between the mass and the radius that depends on the composition of the planet (i.e. rocky or Neptunian) to predict either the mass or radius of an exoplanet for which only one of those values is known (Chen & Kipping 2016).

The measured masses and radii of these planets with their respective errors are shown in Table 1. In order to account for many possible scenarios for the interior structure, 10 000 mass-radius pairs were generated for each planet based on the measured parameters and their errors. For symmetric errors, that is upper and lower errors that differed by less than 20% (i.e. $\left(\left|{\text{Δ}}_{+}-{\text{Δ}}_{-}\right|\right)/\left(\frac{1}{2}\left({\text{Δ}}_{+}+{\text{Δ}}_{-}\right)\right)\le 0.2)$, we applied the same method as Chen & Kipping (2016), who generated random values for masses and radii assuming a Gaussian distribution, with the mean of the errors defined as the standard deviation. The reported masses or radii of all the planets investigated in this study had symmetric errors. Unfortunately, for most of the planets either the radius or the mass was unknown. For these planets, the Forecaster model was used to estimate the planet masses and radii. Since the relation used to predict the missing parameter differed for rocky and Neptunian exoplanets, only mass-radius pairs that implied a rocky composition were used as input parameters for interior modelling. In this study, only those exoplanets that had a probability of more than 5% of having a rocky composition – in other words, where at least 500 of the resulting mass-radius pairs suggested rocky planets – were considered. This narrowed our sample down to a total of 28 exoplanets from the initial 59 HZ rocky planets based on the data from PHL.

The values of planetary eccentricities and semi-major axes were also randomized similarly. The same approach was used in the randomization of these parameters, as described above. In some cases, where only the upper limit of the eccentricity was known, we assumed symmetric errors centred on the measured value, with the upper limit taken as the standard deviation. Only those scenarios that resulted in e ≥ 0 were considered. For each run, all of these parameters were randomly generated in a Monte Carlo approach, which resulted in 10000 distinct combinations of different masses, radii, eccentricities, and semi-major axes for each planet.

2.2 Interior structure modelling

In order to estimate the interior structures of these exoplanets, we used the model described by Dobos et al. (2019), which is a modified version of the model of Barr et al. (2018). In this model planets are assumed to consist of four different layers: an iron core, a rocky mantle, a layer composed of high pressure ice polymorphs (HPP), and a surface ice/water layer. The iron core is assumed to have a uniform density, while the mantle is characterised by compressed Bridgmanite (MgSiO₃). Above pressures of 209 MPa, H₂O is represented as an HPP layer and is assumed to be a mixture of different ice polymorphs from ice II to VII. We used the posterior distributions generated with the Forecaster model where either the planetary mass or radius was unknown. Here we provide a short description of the interior structure model.

For each of the mass-radius pairs, internal compositions that satisfy the following equations were computed: $${\text{Φ}}_{\text{I}}+{\text{Φ}}_{\text{HPP}}+{\text{Φ}}_{\text{rock}}+{\text{Φ}}_{\text{Fe}}=1,$$(1) $${\text{Φ}}_{\text{I}}{\rho }_{\text{I}}+{\text{Φ}}_{\text{HPP}}{\rho }_{\text{HPP}}+{\text{Φ}}_{\text{rock}}{\rho }_{\text{rock}}+{\text{Φ}}_{\text{Fe}}{\rho }_{\text{Fe}}=\rho ,$$(2)

where Φ_I, Φ_ηρρ, Φ_rock, and Φ_Fe are the volume fractions of the surface ice layer, the HPP layer, the rocky mantle, and the iron core, respectively, while ρ_I, ρ_HPP, ρ_rock, and ρ_Fe denote the densities of the same components, respectively. First, the thickness of the ice I layer was estimated. This was achieved by calculating the maximum depth of the surface ice layer (z_I) as $$z_{\text{I}}=P_{\text{I}}/\left({\rho }_{\text{I}}g\right),$$(3)

where P_I = 209MPa is the pressure at which the phase transition of ice I to high-pressure polymorphs occurs, corresponding to a temperature of T = 253 K at the base of the ice I layer. ρ_I = 1000 kg/m³ is the density of the ice layer and $g=G{M}_{\text{pl}}/R_{\text{pl}}^2$ is the surface gravity. Once the volume fraction of the ice I layer was determined, the mean density of the remaining components was calculated. Since no measurements of the sizes of the mantle, iron core, or possible ice layers are available for these planets, the volume fractions of the remaining layers were treated as free parameters. Given the absence of a more thorough understanding of the interior of the modelled exoplanets, uniform probability distribution was assumed for all internal structures in the model. The relative sizes of the HPP layer, rocky mantle, and the iron core were estimated such that their volume fractions satisfied Eqs. (1) and (2). Given the large error bars in the planetary masses and radii (and therefore the densities of the planets), our modelling results in large numbers of possible interior structures with equal probabilities.

2.3 Tidal heating

In order to examine the phase of the ice layer, in other words whether the underlying HPP layer could be in liquid form, we took into account the effect of different heat sources. We calculated the internal heat production of the planets due to heat arising from tidal effects of the host star and due to the decay of radioactive elements. To calculate the planet’s response to tidal forces, a Maxwell viscoelastic rheology was assumed, characterized by a volume-averaged rigidity and viscosity, as described by Dobos et al. (2019). Because of the large uncertainties in the planetary parameters, this approach is more preferable in this present study, since the applied rheology has fewer parameters than more complex rheologies, which yield slightly different behaviors (Henning et al. 2009; Efroimsky 2012). The tidal heat flux was calculated by following the formula of Segatz et al. (1988): $$F_{\text{tidal}}=\frac{21}{8\pi R^2}\text{Im}\left(k_2^{\ast }\right)\frac{R^5{\omega }^5{e}^2}{G}.$$(4)

Here R is the planetary radius, ω = 2π/Ρ is the orbital frequency, e is the eccentricity, G is the gravitational constant, and Im$\left(k_2^{*}\right)$ is the imaginary part of the Love number. This expression is only valid for synchronously rotating bodies with non-zero eccentricities. Assuming a Maxwell viscoelastic rheology, Im$\left(k_2^{*}\right)$ was calculated as $$\text{Im}\left(k_2^{\ast }\right)=\frac{57\eta \omega }{4\rho gR\left[1+{\left(1+\frac{19\mu }{2\rho gR}\right)}^2\frac{{\eta }^2{\omega }^2}{{\mu }^2}\right]},$$(5)

where ρ is the density, μ is the volume-averaged shear modulus, and η is the volume-averaged viscosity. In each layer, the values of both μ and η depend on the state of the given material. As the temperature rises, the fraction of melted material in the layer increases. In response, the viscosity and rigidity of the material decreases. Therefore, both the viscosity and the shear modulus are strongly dependent on temperature. Expressions for the shear moduli and viscosity values of the different layers can be found in Barr et al. (2018).

For the viscoelastic properties of the rock layer, the model described in Dobos & Turner (2015) was used, which is based on previous works of Fischer & Spohn (1990), Moore (2003), and Henning et al. (2009). Here, the tidal response of rock is determined by three threshold temperatures, the solidus temperature (T_s = 1600K), above which the rock becomes partially melted, the breakdown temperature (T_b = 1800K), where the volume fraction of solid rock crystal and melted rock are equal, and the liquidus temperature (T = 2000K), above which the rock is completely melted. This results in four different regimes for the rock layer (below T_s, between T_s and T_b, between T_b and T₁, and above T_l), where the shear modulus and viscosity are expressed differently. A more detailed description of the behavior of the rock layer can be found in Barr et al. (2018). For the HPP and ice I layers, the formulas for viscosity and rigidity also vary based on the state of the layer. If the temperature of the ice layers exceeded the melting point, both μ and η took on the properties of liquid water (μ = 0 , η = 10⁻³ Pa s).

In order to calculate the response of the planet to tidal forces, the properties of the materials in the different layers were taken into account through determining the viscosity and rigidity of each layer and weighting them by their volume fractions (Barr et al. 2018). Similarly to Eq. (2), the averaged viscosity and shear modulus was approximated as $${\text{Φ}}_{\text{I}}{\eta }_{\text{I}}+{\text{Φ}}_{\text{HPP}}{\eta }_{\text{HPP}}+{\text{Φ}}_{\text{rock}}{\eta }_{\text{rock}}\approx \eta ,$$(6) $${\text{Φ}}_{\text{I}}{\mu }_{\text{I}}+{\text{Φ}}_{\text{HPP}}{\mu }_{\text{HPP}}+{\text{Φ}}_{\text{rock}}{\mu }_{\text{rock}}\approx \mu .$$(7)

The planet was then treated as a homogenous body, with the volume-averaged viscosity and rigidity describing its rheolog-ical properties. Therefore, despite using only averaged values for η and μ to determine the tidal heat flux in Eq. (4), in reality the effects of the rock, HPP, and ice layers were all taken into account and incorporated in the volume-averaged viscosity and rigidity. Median values of the volume-averaged rheological parameters are shown in the appendix in Table A.1. By assuming a uniform heat distribution depending only on the radius of the planet, tidal heat flux could be evaluated at any arbitrary distance from the planetary centre. In order to calculate the temperature of the HPP layer, the heat flux generated by tidal forces was first calculated for the surface using the method described above, then scaled to the rock-HPP boundary with a factor of (R/R_rock–HPP)², where R_rock–HPP is the distance of the boundary layer from the centre of the planet.

2.4 Radiogenic heat calculation

In order to evaluate the interior heat flux, the energy released by the decay of radioactive isotopes with long half-lives was also taken into account. We assumed that radioactive elements are only present in the rock mantle and that these elements are the sole source of radiogenic heating. To calculate the radiogenic heat production rate, we used the formula of Hussmann et al. (2010), $$\dot{E}=m_{\text{sil}}{\displaystyle \sum }_{i=1}^4{c}_i{C}_i{H}_i{}_{\exp }\left[{\lambda }_i\left(t-t_{\text{pr}}\right)\right],$$(8)

with i indicating the four long-lived isotopes: ²³⁸U, ²³⁵U, ²³²Th, and ⁴⁰K. The mass of the rock mantle, m_sil = M_pl × M_rock, where M_pl is the mass of the planet and M_rock is the mass fraction of the rock mantle, determined with the interior model (see Sect. 2.2). The values used for the composition, c, and heat release, H, of the different elements are shown in Table 2 (Audi et al. 1997). Since we did not have information about the abundances (C) of these isotopes in the investigated systems, these values are based on different types of meteorites found on Earth (Lodders & Fegley 1998). We considered two extreme cases to estimate the elemental abundances: CI and LL chondrite abundances. As a lower limit we chose CI chondrite abundances, since this produces the lowest energy rates from the sample. For similar reasons, LL chondrite abundances were chosen as the higher limits. The elemental abundances for each radioactive isotope were then approximated as C = fC_CI + (1 − f /100)C_LL, where f was an integer between 0 and 100. The value for f was randomized with uniform probability for each internal structure. In Eq. (8), λ_i is the decay constant, and depends on the half-lives, τ_i, as λ_i = in 0.5/τ_i. The values of the half-lives and the elementa abundances are shown for all four isotopes in Table 2. Since our calculations focused on the current heat production of the exoplanets, here the parameter, t, denotes the age of the system, while t_pr = 4.6 × 10⁹ Gyr is the age of the Solar System. In the case of GJ 273, K2-72, and Kepler-1649, where the ages of the stars were unknown, t = t_pr was used. The radiogenic heat flux, F_rad, was evaluated at the boundary of the rock mantle and the HPP layer using $$F_{\text{rad}}=\frac{\dot{E}}{4\pi R_{\text{rock–HPP}}}.$$(9)

2.5 Internal heat dissipation

We assumed that the heat produced by both tidal interactions and the decay of radioactive isotopes was transported out of the system by solid-state convection. The temperatures at the top and bottom of each convective layer are constant in space due to the 1D nature of our model, but vary as a function of the heat flow across each layer. This accounts for heat transport across the compositional or phase boundary layers by conduction (e.g. Mueller & McKinnon, 1988). The convective flux was calculated by using the formula of Barr (2008): $$F_{\text{conv}}=0.53k{\left(\frac{\rho g\alpha }{\kappa }\right)}^{1/3}{\left[\frac{Q}{R_{\text{G}}{T}_{\text{hpp}}^2}\right]}^{-4/3}{\eta }_{\text{hpp}}{\left(T_{\text{hpp}}\right)}^{-1/3},$$(10)

where k is the thermal conductivity, ρ is the density of the HPP layer, g is the local gravitational acceleration, α is the thermal expansion coefficient, κ is the thermal diffusivity, Q is the activation energy for ice, R_G = 8.314 J mol⁻¹ K is the gas constant, and T_hpp is the temperature of the HPP layer in K. The viscosity of the HPP layer, η_hpp, can be expressed as η_hpp = c_rh exp [Q/(RGT_hpp)] by following the flow law of Showman & Han (2005), where c_rh = 10¹⁸/ exp [Q/(R_GT_melt,hpp)] is the rheological constant for the high-pressure ice polymorph, and T_melt,hpp is the melting temperature of the high-pressure ice polymorph, which was assumed to be 280 K in this work. This is a reasonable assumption for the whole layer that spans from the bottom of the ice I layer at 209 MPa to pressures around 2.5 GPa, where the Ice VII phase dominates and above which the rheology of ice is not well constrained. In reality, the melting temperature of the HPP is the lowest at the top of the HPP layer (253 K) and the highest at the lower boundary. The values of the physical parameters of the HPP layer are shown in Table 3.

The equilibrium temperature of the HPP layer, Thpp, was calculated by equating the convective flux with the internal heat flux at the rock-HPP boundary: $$F_{\text{conv}}=F_{\text{int}}=F_{\text{tid}}+F_{\text{rad}},$$(11)

where both tidal and radiogenic heat fluxes represent values scaled by the area of the rock-HPP boundary surface (see Sect. 2.3.). We assumed that H₂O was present in liquid form underneath the surface layer, when the equilibrium temperature of the HPP layer was above the melting point of the high-pressure ice polymorph.

2.6 Summary of the analysis

Here we present a brief summary of how we applied the methods detailed in the previous sections. We modelled the interiors of 28 HZ exoplanets based on their measured parameters. Planetary masses and radii were randomly generated with the Forecaster model using an empirical relation, if either one was unknown. Since a large number of interior structures were possible for a given density, the generated mass-radius pairs corresponded to tens to hundreds of thousands of distinct possible interiors per exoplanet.

We calculated the total internal heat flux at the rock-HPP boundary for each internal composition. We assumed two different heat sources: tidal heating and radiogenic heating. For the latter, elemental abundances were approximated using values measured in CI and LL chondrite meteorites. We note that other heat sources, such as secular cooling of the core or greenhouse warming above the surface may be additional contributors to the heating of the interior. Our calculations therefore serve as lower estimates for the total internal heat flux. We calculated the temperature of the HPP layer by assuming that the internal heat flux was transported out of the system via solid-state convection. If the equilibrium temperature of the HPP layer was above the melting point of the HPP ice, we concluded that the HPP layer was in liquid form. This was calculated for all the possible interior structures of each planet.

3 Results

3.1 H₂O content of the planets

We estimated the H₂O content of all 28 planets based on their modelled interiors. For a given internal structure, the H₂O mass fraction is defined as the sum of the mass fractions of the surface ice/water layer and the underlying HPP layer. Since a large number of different internal structures were feasible for a given planet, our calculations resulted in a wide range of H₂O mass fractions for each planet.

In Fig. 1 the distribution of the resulting internal structures in different H₂O mass fraction ranges for each planet is shown. Each box represents a 3% mass fraction interval in the top panel of the figure and 6% in the bottom panel, respectively. The colour and the numbers indicate the fraction of all modelled interiors that resulted in internal structures within a certain H₂O mass fraction range. If a large number of interiors had high H₂O mass fractions for a given planet, that appears as a shift in the values and lighter colours in our figure towards the higher mass fraction ranges. If the input parameters - the mass and radius of the planet - had low uncertainties, the resulting H₂O mass fractions concentrate in a relatively narrow range, with the H₂O mass fractions showing a peak in the distribution at the largest probabilities. For planets with less constrained parameters, the H₂O mass fractions spanned a broader range of feasible values (see for example K2-3 d, K2-72 e or Kepler-1229 b).

Our results show that all of these planets could have large enough H₂O mass fractions to have global ice/water surfaces. Planets with surface temperatures above the melting point of ice and large water mass fractions (e.g. Kepler-296 e) may have global oceans regardless of the presence of an atmosphere. The equilibrium temperatures (T_eq) of the investigated systems are shown in Table 1. It is important to note that the method for estimating T_eq may differ for the individual planets in the assumed heat distribution and Bond albedos. For more details, the reader is referred to the papers referenced in Table 1. For those planets where the surface temperature falls below the melting point of ice, liquid water may not be present on the surface in the absence of an atmosphere, and large water mass fractions may imply global ice layers. If the internal heat flux in these bodies is sufficiently high to cause melting in the HPP layer, these worlds may harbour underground liquid water reservoirs. Interactions between the rock and water layer may then result in a suitable environment for life to develop, which makes these planets interesting targets for habitability studies. Since the investigated planets all lie within the HZ, high values for H₂O mass fractions could also imply global oceans on the surface if the planets have atmospheres similar to that of the Earth.

For some planets, like GJ 1061 c and d, the thickness of the ice layer could not be properly determined, since our model resulted in a similar number of cases in a wide range of water mass fractions. This arose from the fact that the planetary parameters of these bodies were not known to sufficiently high precision. In order to refine our results for these planets, more precise measurements of the planetary masses and radii are required.

For planets with well-constrained planetary masses and radii, such as the TRAPPIST-1 system, the interior structures were constrained to a higher degree, which results in more precise H₂O-content estimations. All of the modelled TRAPPIST-1 planets are likely to have extended H₂O layers. Since TRAPPIST-1 e, f, g, and potentially planet d as well are in the HZ of their star, determining the state of the H₂O layer is necessary to assess the habitability of these planets.

In some cases the modelled exoplanets had extremely large H₂O mass fractions, which may imply that those planets fall under the category of ocean worlds. Although our model did not account for planetary atmospheres, HZ exoplanets with such high water mass fractions are likely to have surface oceans, since water vapour will provide non-negligible atmospheric pressure. Our results suggest that Kepler-62 f, Kepler-452 b and Kepler-442 b could be members of the ocean world class (bottom panel of Fig. 1). In the case of Kepler-62 e and K2-288 B b, most of the generated mass-radius pairs result in mini-Neptunian compositions. However, if indeed superterrans, our models show that these planets may also be part of the ocean world planetary group. For each planet, the probabilities of having rocky compositions (here referred to as terran probability) are shown in the seventh column of Table 4.

Fig. 1 Distributions of the modelled H2O contents of the investigated exoplanets. Here the H2O mass fractions are divided into 3% or 6% groups in the top and bottom panels, respectively. The colour indicates the fraction of internal structures within the given percentage groups of the water mass fraction, with the values written on the figure. Probabilities of less than 1% are not shown. Planets with water mass fractions over 30% are shown in the bottom panel. We note the different scales in the two panels.

3.2 Effects of tidal and radiogenic heating

The separate effects of radiogenic and tidal heating at the rock-HPP ice boundary are shown in Fig. 2 and in Table 4. The median values of radiogenic and tidal heat fluxes are shown in Table 4 along with values corresponding to the first and third quartiles. Boxplots for tidal and radiogenic heating are shown in Fig. 2. The top panel depicts planets with strong tidal heating, while the bottom panel shows exoplanets where radiogenic heat flux is comparable with or exceeds the flux from tidal heating. The maximum values correspond to Q3 + 1.5 × IQR, while the minim are defined as Q1 − 1.5 × IQR, where IQR = Q3 − Q1 is the interquartile range. In the case of radiogenic heating, elemental abundances were in every case a mix of CI and LL chondritic abundances. Tidal heating, if it was present in the planet, was almost always at least an order of magnitude higher than radiogenic heating and was the main source of internal heat for close-in exoplanets around late-type stars, as seen in the top panel of Fig. 2. Tidal heat flux on these planets, like GJ 1061 c and d, was close to 2 W m⁻², similar to the heat flux observed on Io (Veeder et al. 2004). In some planets, as in the case of Kepler-296 e, the average heat flux from the two different sources is comparable.

On average, the internal heat flux provided by radiogenic heating was significantly lower compared to the tidal heat flux in close-in exoplanets. Melting of the HPP ice layer, however, still took place even in the absence of tidal interactions on many of the investigated exoplanets, as seen in Table 4. In the case of some exoplanets with radiogenic heating as the main internal heat source, the assumed abundances for radioactive elements may be the decisive factor for melting to occur. More data on radioactive elemental abundances in exoplanetary systems may improve our understanding of the interior heating in these bodies.

Planets with considerable tidal heating (shown in the top panel of Fig. 2) were further investigated, in order to have a better insight into the possible tidal heating rates. For these planets the distribution of the tidal heat flux at the rock-HPP boundary is shown in Fig. 3. Tidal heat flux values shown here were averaged for all mass-radius pairs in the case of every planet, which gave a total of maximum 10 000 different tidal heat flux values. Structures where tidal heating had no contribution to the internal heat budget were possible for all exoplanets, though to different degrees for each body (represented by the height of the peaks at zero W m⁻²). The shape of the distribution, as well as the probable heat flux values, also differed for the investigated planets. It is clearly seen from the histograms that the maximum values for the tidal heating rates of many of these planets (e.g. GJ 1061 c and d) were not representative values, but were in fact outliers in the distribution.

Fig. 2 Contribution of the different heat sources to the internal heat flux at the rock-HPP boundary for each planet. The red boxes indicate tidal heating, while blue denotes radiogenic heat flux. Median values are represented as vertical lines inside the boxes. The upper and lower boundaries of the boxes correspond to the third and first quartiles (Q3 and Q1), respectively, while error bars show the maximum and minimum values. White stars indicate the average tidal heat flux values for every planet. We note the different scales in the two panels.

3.3 Melting caused by internal heating

The phase of H₂O in the HPP layer was further investigated. If the combined internal heat flux from tidal heating and the decay of radioactive nuclei for a given internal structure was high enough to heat up the HPP layer above the melting point, it resulted in a water layer below the surface ice. Here, the melting probability is the fraction of all modelled structures that resulted in T_HPP >T_melt,HPP. Melting probabilities for each planet are shown in Table 4. The higher the melting probability for a certain planet, the more likely that the HPP layer is melted.

The effects of tidal and radiogenic heating on the melting of the HPP layer were studied in further detail. For close-in HZ exoplanets, for example GJ 1061 c and TRAPPIST-1 d, tidal heating was efficient enough to melt the underlying ice layer in most of the cases. This caused a higher overall fraction of modelled structures to experience melting than in the case of exoplanets farther out from their host stars, where radiogenic heating was the main contributor to the internal heat budget. For planets in the HZs of more massive stars, tidal dissipation is not a substantial source of energy, since for these planets the probability of being captured into spin-orbit resonance in small.

Since the radiogenic heat production rate depends on the thickness of the rock mantle, different internal structures resulted in different internal heat fluxes. Differences in the internal structures therefore affected the resulting radiogenic heat flux. Further uncertainties in the melting probabilities originated from the unknown abundances of radioactive elements in the investigated exoplanetary systems. In this study these abundances were approximated based on a mixture of elemental abundances seen in Solar System meteorites, with the two extremes being the purely LL chondritic abundances (providing the maximum possible heat flux for exoplanets dominated by radiogenic heating) and the purely CI chondritic abundances (associated with the lowest heating rates, see Sect. 2.4).

Fig. 3 Histograms showing the distribution of the tidal heat flux values, Ftidal, calculated at the rock-HPP boundary for planets with strong tidal heating, with N denoting number density.

3.4 Overall likelihood of containing liquid water

In order to properly assess the potential of a water reservoir under the surface of these planets, it is important to take into account other factors, such as the likelihood of the planet having a rocky composition. The values for terran probability, in other words the probability of a planet having a rocky composition, are represented by the fraction of generated planetary masses below 2.04 M_⊕ for each planet. This boundary mass is based on the transition point between rocky and Neptune-like planets in the Forecaster model. On the other hand, in a recent study Luque & Pallé (2022) argue that for planets around M dwarf stars, the density of the exoplanet is the key parameter that separates rocky planets, water worlds, and mini-Neptunes from each other, not the planetary mass or radius. However, for the sake of consistency, in this present study the boundary mass based on the transition point in the Forecaster code was used for classification in all cases. The only exception was TOI-1452 b, where both the planetary mass and radius with symmetric errors were known. In this case, even though the mass of the planet was over the 2.04 M_⊕ transition point, the density suggested a terrestrial composition in all cases. Since neither the mass nor the radius had to be approximated using the Forecaster code, this did not result in any inconsistency. It is important to note that while some planets, for example Kepler-442 b, may have an extended H₂O-layer that is likely to melt due to internal heating processes, the planet itself has a low probability of being a rocky planet. It is evident that the best targets for habitability studies are those exoplanets that are highly likely to have rocky composition, have a high chance of having an H₂O-layer, and experience moderate internal heat flux.

To identify the most likely candidates for having a liquid water layer, we defined the ocean probability parameter (P_ocean) as $$P_{\text{ocean}}=P_{\text{terran}}\ast P_{\text{melt}},$$(12)

where P_terran is the terran probability and P_melt is the melting probability defined in Sect. 3.3. The P_ocean parameters are shown in Table 4 for each of the investigated planets. Our calculations show that Ross 128 b, Teegarden’s b, Proxima Centauri b, TRAPPIST-1 d, GJ 1061 c and d, Teegarden’s c, and TRAPPIST-1 e have the largest values for P_ocean, and hence are the most likely to have extended internal water reservoirs.

Figure 4 shows P_ocean for the investigated systems, along with the HZ for stars with effective temperatures above 2750 K. The HZ is shown for illustrative purposes only. The boundaries of the HZ were calculated based on the stellar parameters (mass, temperature, and luminosity) of 2.5 Gyr old stars from the MESA Isochrones and Stellar Tracks (MIST) online tool (Dotter 2016; Choi et al. 2016; Paxton et al. 2018). A second-order polynomial was fitted to the data using the method of Kopparapu et al. (2014) for planets with masses of 1 Μ_⊕, which provided the boundary of the HZ for stars above 2600 K. The HZ of TRAPPIST-1 is not shown in the figure because the temperature of the star is below this value, and this method could not be used to estimate its HZ.

4 Discussion and conclusions

Our research provides an insight into the conditions of the interiors of 28 HZ rocky exoplanets from the viewpoint of habitability. We modelled the interior structure of each planet and determined the possible H₂O content of these bodies. Internal heat fluxes from the combined effect of tidal and radiogenic heating were calculated for all exoplanets. We estimated the temperature of the HPP ice layer based on the internal heat flux and calculated the probability of the HPP layer melting. Our calculations suggest that planets with the highest probabilities of containing liquid water reservoirs are Teegarden’s b and c, GJ 1061 c and d, TRAPPIST-1 d and e, Ross 128 b, and Proxima Centauri b. These exoplanets are likely to be rocky planets while having an extended water layer, and hence are of outstanding importance to habitability studies.

An interesting remark is that all these planets orbit M dwarf stars. This, however, does not mean that exoplanets of M dwarfs are more likely to have subsurface oceans in general, but rather it implies that there is some kind of an underlying bias in our sample. A possible explanation is that our current exoplanet detection methods are more sensitive to finding planets with lower masses or of smaller sizes orbiting closer to their host stars than rocky worlds of similar sizes further out in their systems. Hence, HZ exoplanets with smaller sizes or masses are more likely to be detected around M dwarf stars. Since P_terran is determined based on planetary masses, low-mass bodies are expected to have larger probabilities of having rocky compositions than planets with higher measured sizes or masses. This could explain why planets in our sample that reside in the HZ of earlier-type (F,G,K) stars have smaller P_ocean – a parameter that is directly proportional with P_terran – than exoplanets orbiting M dwarf stars.

There were a number of exoplanets in our sample (such as GJ 1061 c and d) where the internal heating processes proved to be strong enough to induce melting in more than 98% of all modelled interiors. These planets are of particular interest, since melting is almost guaranteed to take place. With the currently available mass and radius measurements, however, it is not possible to properly assess the probability of a thick H₂O layer being present on these planets. Although the results of our modelling permitted the presence of an extended ice layer, the H₂O mass fraction of the planet could not be adequately constrained. To have a greater insight into the structures of these potentially habitable bodies, more precise observations of their masses and radii are required.

In order to properly estimate the radiogenic heat flux, a better understanding of the average elemental composition in the investigated systems is required. A more realistic approximation of the radioactive isotope abundances may improve our knowledge of the internal heat production in these exoplanets. In accordance with previous studies on this subject, which showed that solar analogues may be more abundant in radioactive elements than our Sun, we assume that our calculations do not overestimate the radiogenic heat flux in the investigated systems.

It is worth noting that tidal heating is strongly dependent on the eccentricity of the planetary orbit, which changes during the evolution of the planet. As a result of tidal forces, the orbit of close-in exoplanets may eventually become circular (e = 0). Tidal heating may therefore completely diminish without any external perturbation that would help maintain non-zero eccentricities. However, e > 0 may prevail due to the gravitational effects of other planets in the system. For example, in the TRAPPIST-1 system, mean motion resonances prevented the circulization of planetary orbits. The eccentricities of the TRAPPIST-1 planets are still above zero, despite the old age of the system.

Since the atmosphere contributes little to the mean density of rocky planets compared to the mantle or the iron core (Noack et al. 2017), for our approach of modelling the interiors of exo-planets the presence of an atmosphere was not assumed. The knowledge of whether the planet has an extended atmosphere or not, however, could help in the interpretation of our results. One such planet is Kepler-62 f, where an extended, low-density atmosphere could explain the lower mean density of the planet, which in our modelling was interpreted as extremely high water content. Furthermore, liquid oceans on the surface of rocky planets with low water mass fractions can be present if the planets have atmospheres. On the other hand, planets with high water mass fractions and equilibrium surface temperatures above the melting point of ice may have global oceans on their surfaces even without an extended atmosphere and regardless of internal heating. TRAPPIST-1 d and Kepler-296 e may be such bodies, where the evaporation of a possible surface ocean itself may be sufficient to produce a low-pressure water vapour atmosphere. However, planets without a surface ice layer, and with internal heat flux high enough to induce a runaway greenhouse state, could end up losing their water (Barnes et al. 2013). The effect of tidal heating on the habitability of exoplanets and its possible impact on the boundaries of the HZ is further explored in Barnes et al. (2009).

Modelling the interior of these rocky bodies gives us a more accurate picture of the structure and internal conditions of HZ exoplanets. Determining the H₂O content and the internal heating of these planets allows us to estimate the probability of a large water reservoir and opens a new window for habitability studies. Subsurface oceans may be excellent sites for the development of life, although this may be difficult to observe with our current technology.

According to Tian & Ida (2015), planets with Earth-like water content around M dwarf stars may be rare because of water loss due to the high XUV radiation of the star. However, they found that ocean planets with a water mass fraction higher than 10% are able to retain their water, a fact that coincides with our results on the water-rich interiors of the TRAPPIST-1 planets. Recent studies argue that the origin of high water content in planets around low-mass stars can be explained by migration models (Unterborn et al. 2018; Miguel et al. 2020). Unterborn et al. (2018) found that the TRAPPIST-1 planets may have formed outside of the snow line and then migrated inward. On the other hand, water transport through asteroid bombardment cannot be ruled out as well. Bombardment in compact systems, which could be responsible for a smaller fraction of the total water content observed in these planets, was found to be possible through various methods, for example assuming a giant planet within the system or a perturbing body outside of the stellar system (Ðošovic et al. 2020; Clement et al. 2022).

Brugger et al. (2016) modelled the structure of Proxima Cen-tauri b assuming five different layers – a metallic core, a lower and an upper mantle, a HPP layer, and a surface H₂O layer. They concluded that even internal structures with water mass fractions of 50% are possible for Proxima Centauri b. Herath et al. (2021) investigated the interior structures of Ross 128 b and Proxima Centauri b and – in agreement with Brugger et al. (2016) for the latter and with our results for both exoplanets – showed that high volatile content in these planets cannot be excluded. They investigated scenarios with different water mass fractions in the case of both planets, and concluded that a liquid water layer can be present over a thick ice shell for structures with water mass fractions larger than 5%, assuming a planetary atmosphere. Our results also suggest that there is a high chance that Proxima Centauri b may have a liquid ocean layer, although, contrary to Brugger et al. (2016), in our study this liquid water layer is beneath the ice shell. Nevertheless, in the presence of an atmosphere, our results are consistent with a surface ocean layer. Furthermore, Brugger et al. (2016) also find that high water content can hinder the generation of a magnetic field in these planets. Notably, they conclude that water mass fractions larger than 10% make the presence of liquid outer cores unlikely.

Interior modelling for TOI-1452 b by Cadieux et al. (2022) suggests that the planet may be a water world with a high water mass fraction of ${22}_{-13}^{+21}\%$. In our study, however, scenarios with a significantly lower water mass fraction are favoured, although high values like those presented by Cadieux et al. (2022) cannot be excluded either.

In agreement with our study, Jackson et al. (2008a) found that tidal heating may be beneficial for the subsurface habitability of icy bodies by allowing the existence of an ocean under the surface ice layer.

Besides modelling the interior of HZ exoplanets, our approach is also suitable for the characterization of rocky exo-planets farther outside the HZs of their stars. However, the only

"cold" rocky exoplanets known to date are those discovered by gravitational microlensing. Once such exoplanets with the potential of further characterization are discovered, our method could be applied to gain more insight into the internal conditions of those bodies.

Fig. 4 Internal ocean probabilities of the investigated exoplanets. The colour bar depicts the probability of each planet having a liquid water reservoir in its interior. Red colours imply a high probability, while planets with dark blue circles are less likely to have internal oceans. The HZ, shown as a green band only for illustrative purposes, was calculated based on stellar parameters from theoretical stellar evolution models.

Acknowledgements

The authors thank the help of Zoltán Dencs, Michaela Walterová, Eike Guenther and René Heller for discussing the results of the manuscript. Á.B. and V.D. have been supported by the Hungarian National Research, Development, and Innovation Office (NKFIH) grant K-131508. The COFUND project oLife has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement no. 847675. This research has made use of the Planet Habitability Laboratory website hosted by the University of Puerto Rico at Arecibo.

Appendix A Volume-averaged rheological parameters

Median values for the volume-averaged shear modulus, μ, and viscosity, η, are shown, along with the values of the first and third quartiles in the lower and upper indices, respectively.

References

All Tables

All Figures

Fig. 1 Distributions of the modelled H2O contents of the investigated exoplanets. Here the H2O mass fractions are divided into 3% or 6% groups in the top and bottom panels, respectively. The colour indicates the fraction of internal structures within the given percentage groups of the water mass fraction, with the values written on the figure. Probabilities of less than 1% are not shown. Planets with water mass fractions over 30% are shown in the bottom panel. We note the different scales in the two panels.

In the text

Fig. 2 Contribution of the different heat sources to the internal heat flux at the rock-HPP boundary for each planet. The red boxes indicate tidal heating, while blue denotes radiogenic heat flux. Median values are represented as vertical lines inside the boxes. The upper and lower boundaries of the boxes correspond to the third and first quartiles (Q3 and Q1), respectively, while error bars show the maximum and minimum values. White stars indicate the average tidal heat flux values for every planet. We note the different scales in the two panels.

In the text

	Fig. 3 Histograms showing the distribution of the tidal heat flux values, Ftidal, calculated at the rock-HPP boundary for planets with strong tidal heating, with N denoting number density.
In the text

Fig. 4 Internal ocean probabilities of the investigated exoplanets. The colour bar depicts the probability of each planet having a liquid water reservoir in its interior. Red colours imply a high probability, while planets with dark blue circles are less likely to have internal oceans. The HZ, shown as a green band only for illustrative purposes, was calculated based on stellar parameters from theoretical stellar evolution models.

In the text