Potential for Liquid Water Biochemistry Deep under the Surfaces of the Moon, Mars, and beyond

The habitable zone (HZ), in its classical form, delineates the region around the host star where liquid water could exist on the surface of a rocky planet (Kasting et al. 1993). The concept of the HZ has been extended in myriad directions by including new greenhouse gases, taking surface geology into account, expanding it to the galactic level, and many more (Ramirez 2018), thereby making it widespread in astrobiology, especially from a practical standpoint in the search for biosignatures. It is worth emphasizing, however that the canonical HZ (Kasting et al. 1993) definition pertains to remote detection of habitability. Worlds outside the HZ may very well be habitable, as illustrated by the classic examples of Europa, Enceladus, and Titan in our solar system. Another notable category of potentially habitable worlds is composed of free-floating planets and moons, which might either host liquid water on the surface (Stevenson 1999; Lingam & Loeb 2020) or comprise oceans beneath icy envelopes (Abbot & Switzer 2011).

In a similar vein, one may investigate the habitability of rocky environments underneath the surface. On Earth, subsurface rock-based (endolithic) environments are known to host rich and thriving biospheres (Gold 1992; Edwards et al. 2012), and the total biomass in these settings is predicted to be ∼10% of Earth's overall value (Whitman et al. 1998; Bar-On et al. 2018). Hence, it is not surprising that a number of studies have investigated Martian subsurface habitability (Boston et al. 1992; Cockell 2014; Michalski et al. 2018; Sholes et al. 2019; Carrier et al. 2020). However, generic habitability studies of deep terrestrial biospheres are lacking, with the exception of McMahon et al. (2013).

Our work differs from the aforementioned publication in the following respects: (i) we explicitly account for the possibility of free-floating rocky objects with only internal heating, (ii) we take other constraints (e.g., pressure and energy) into consideration, and (iii) we span a larger parameter space of physical variables. We calculate the extent of the habitable subsurface region in Section 2, tackle other properties of deep biospheres in Section 3, and summarize our findings in Section 4.

The minimum temperature for liquid water is roughly between 250-270 K when the pressure is ≲1 GPa (Choukroun & Grasset 2007). The cryophile Planococcus halocryophilus grows at 258 K (Merino et al. 2019), which falls in this regime. Both experiments and theory indicate that lower temperatures would result in the vitrification of cells. The maximum temperature for liquid water's existence is sensitive to total pressure and exceeds 373 K when the pressure is > 1 bar. However, at sufficiently high temperatures, a number of inexorable biochemical and biophysical constraints come into play (Clarke 2014; McKay 2014). The hyperthermophile Methanopyrus kandleri survives at 395 K (Takai et al. 2008), although the limits for organic life might extend up to ∼450 K (Bains et al. 2015). Based on the preceding considerations, T_min = 260 K and T_max = 400 K demarcate our thermal range.

In its most general form, the heat-transfer equation is expressible as (Melosh 2011, Equation 4.13):

$$\begin{eqnarray}&&\displaystyle \frac{\partial T}{\partial t}={\rm{\nabla }}\cdot \left(\kappa {\rm{\nabla }}T\right)+\displaystyle \frac{{ \mathcal H }}{{C}_{P}},\end{eqnarray} \tag{ 1 }$$

where T is the temperature, whereas κ, C_P and ${ \mathcal H }$ denote the thermal diffusivity, specific heat capacity at constant pressure, and the rate of heating per unit mass. As our goal is to develop a simple model, we will follow the approach delineated in Lingam & Loeb (2019), which generalizes the slab models described in Abbot & Switzer (2011) and McMahon et al. (2013).

To begin with, we tackle the steady-state solution; although ${ \mathcal H }$ is time-dependent, it is subject to significant change only over long timescales of several Gyr (Turcotte & Schubert 2002, Section 6.24). We suppose that Fourier's law is valid, thereby yielding

$$\begin{eqnarray}&&{ \mathcal Q }+k\displaystyle \frac{{dT}}{{dr}}=0,\end{eqnarray} \tag{ 2 }$$

where ${ \mathcal Q }\equiv { \mathcal Q }(r)$ is the heat flux at r and $k\equiv k(T)$ denotes the thermal conductivity of the rock(s). The former quantity is modeled by using

$$\begin{eqnarray}&&{ \mathcal Q }=\displaystyle \frac{Q}{4\pi {r}^{2}}\times \displaystyle \frac{4\pi {r}^{3}/3}{4\pi {R}^{3}/3}=\displaystyle \frac{Q}{4\pi {R}^{2}}\displaystyle \frac{r}{R},\end{eqnarray} \tag{ 3 }$$

where R is the radius of the world, while Q embodies the total internal heat flow. Therefore, the above ansatz posits that the heat flux is the ratio of the heat flow at this location (Q_E) and the area encompassed (4π r²). The former variable is determined by assuming that the heating is uniformly distributed in the interior, amounting to ${Q}_{E}\approx Q\times \left(\tfrac{4\pi }{3}{r}^{3}/\tfrac{4\pi }{3}{R}^{3}\right)$. Furthermore, we specify Q following the standard prescription (e.g., Valencia & O'Connell 2009; Lingam & Loeb 2021):

$$\begin{eqnarray}&&Q={\rm{\Gamma }}\,{Q}_{\oplus }\left(\displaystyle \frac{M}{{M}_{\oplus }}\right),\end{eqnarray} \tag{ 4 }$$

where ${Q}_{\oplus }\approx 4.4\times {10}^{13}$ W is the Earth's heat flow (The KamLAND Collaboration 2011), Γ is the abundance of long-lived radionuclides per unit mass measured relative to Earth, and M is the object's mass, which is related to R via the scaling $M\propto {R}^{3.7}$ (Zeng et al. 2016). Note that Equations (2) and (3) ensure that the geothermal flux condition is satisfied at the surface, and dT/dr = 0 at r = 0.

The remaining component is k(T). Evidently, this quantity is subject to much variability as it depends upon the characteristics of the rocks (e.g., sedimentary, metamorphic). For instance, the thermal conductivity of sandstone is ∼2 times higher than basalt and granite at 273 K (Hartlieb et al. 2016). This issue is compounded by the absence of empirical studies of heat conductivity at very low temperatures, which has necessitated reliance on heuristic models. Bearing these caveats in mind, we will consider a rocky object with a Mars-like composition, whose thermal conductivity at the surface is ≲2 times smaller than that of Earth (cf. Whittington et al. 2009; Parro et al. 2017). For this model, we treat the crust as being dominated by basalt and employ the relationship (Clifford et al. 2010, their Equation 1):

$$\begin{eqnarray}&&k(T)\approx { \mathcal A }+\displaystyle \frac{{ \mathcal B }}{T},\end{eqnarray} \tag{ 5 }$$

where ${ \mathcal A }=0.4685$ W m⁻¹ K⁻¹ and ${ \mathcal B }=488.19$ W m⁻¹. After integrating Equation (2) and drawing upon the preceding formulae, we arrive at

$$\begin{eqnarray}&&\displaystyle \frac{Q\left({R}^{2}-{r}^{2}\right)}{8\pi {R}^{3}}={ \mathcal A }\left(T-{T}_{s}\right)+{ \mathcal B }\mathrm{ln}\left(\displaystyle \frac{T}{{T}_{s}}\right),\end{eqnarray} \tag{ 6 }$$

where T_s is the average surface temperature.³ We define $F(T)\equiv { \mathcal A }\left(T-{T}_{s}\right)+{ \mathcal B }\mathrm{ln}(T/{T}_{s})$, which plays a prominent role in the subsequent analysis.

There are two distinct cases that we must investigate in order to determine the spatial extent H of the putative subsurface biosphere. In the first scenario, we have T_s > T_min so that the conditions are warm enough in principle for surface-based life to survive. In this case, as explained in the Appendix, we find that

$$\begin{eqnarray}&&\displaystyle \frac{H}{R}=1-\sqrt{1-\displaystyle \frac{8\pi R}{Q}F({T}_{\max })},\end{eqnarray} \tag{ 7 }$$

where F(T) was introduced below Equation (6). For $H\ll R$, simplifying Equation (7) further leads to

$$\begin{eqnarray}&&H\approx \displaystyle \frac{4\pi {R}^{2}F({T}_{\max })}{Q}.\end{eqnarray} \tag{ 8 }$$

In the second case, we consider T_s < T_min, implying that life near the surface is not feasible. In this setting, H is computed in the Appendix and found to be

$$\begin{eqnarray}&&\displaystyle \frac{H}{R}=\sqrt{1-\displaystyle \frac{8\pi R}{Q}F({T}_{\min })}-\sqrt{1-\displaystyle \frac{8\pi R}{Q}F({T}_{\max })}.\end{eqnarray} \tag{ 9 }$$

If the lower and upper limits of this biosphere, delineated in the Appendix, are much smaller than R, the above equation simplifies to

$$\begin{eqnarray}&&H\approx \displaystyle \frac{4\pi {R}^{2}}{Q}\left[F({T}_{\max })-F({T}_{\min })\right].\end{eqnarray} \tag{ 10 }$$

For the second case, we note that free-floating rocky objects would fall under this category. The surface temperature is purely set by internal heat and is determined via the Stefan-Boltzmann law, which yields

$$\begin{eqnarray}&&{T}_{s}\approx 35\,{\rm{K}}\,{{\rm{\Gamma }}}^{1/4}{\left(\displaystyle \frac{R}{{R}_{\oplus }}\right)}^{0.43}.\end{eqnarray} \tag{ 11 }$$

From inspecting Equations (7) and (9), we notice that H is real-valued only if the following criterion holds true:

$$\begin{eqnarray}&&\displaystyle \frac{Q}{8\pi R}\gt F({T}_{\max }).\end{eqnarray} \tag{ 12 }$$

In what follows, we will implicitly presume that this condition is fulfilled. Broadly speaking, this translates to $R\gtrsim 0.1\,{R}_{\oplus }$, with the caveat that the exact magnitude is sensitive to T_s.

From the preceding discussion, evidently H is a complex function of Γ, R, and T_s, except for free-floating rocky objects (see the Appendix). In Figure 1, we have plotted H as a function of T_s and R while fixing Γ = 1. Note that this figure corresponds to rocky objects where the surface temperature is set by either stellar or tidal heating and not by internal (i.e., radiogenic and primordial) heating. There are two broad trends that can be inferred. First, H exhibits a fairly strong dependence on R, which becomes evident from Equations (8) and (10) as well. Second, H is nearly independent of T_s for T_s < T_min, viz., H is largely unaffected by the properties of the region between T = T_s and T = T_min.

Zoom In Zoom Out Reset image size

As we chose a Mars-like composition, we compare our results against prior studies of the Martian subsurface biosphere. From the Appendix, we ascertain that the habitable subsurface region spans ${H}_{1}\approx 4.2$ km and H₂ ≈ 13.6 km beneath the surface, thereby corresponding to H = 9.4 km. Note that H₁ and H₂ signify the locations of the top (T = T_min) and bottom (T = T_max) layers of the subsurface habitable region, as measured from the object's surface. This result is consistent with the slab model of McMahon et al. (2013), where it was determined that the biosphere might span H₁ ≈ 4.5 to H₂ ≈ 14.0 km, thus yielding H = 9.5 km. In a more detailed analysis, Clifford et al. (2010) estimated H ≈ 2.3−4.7 km at the equator and H ≈ 6.5–12.5 km at the poles. If one takes the arithmetic mean of these two extremes, we end up with H ≈ 4.4−8.6 km; the upper bound is therefore not far removed from our model's predictions.

Next, we turn our attention to free-floating objects where H is only a function of R and Γ as outlined in the Appendix. We have plotted the ensuing behavior of H in Figure 2 for different choices of Γ. Aside from the expected sensitivity to R, we see that Γ plays a fairly important role in regulating the spatial extent of the habitable region. An interesting point worth mentioning is that higher Γ augurs well for surface habitability in some respects (Lingam & Loeb 2020), but it suppresses the extent of the subsurface habitable zone and other characteristics encountered hereafter because of the inverse dependence on Γ in Equations (8) and (10), after using the scaling for Q prescribed in Equation (4).

Zoom In Zoom Out Reset image size

In this Section, we sketch other constraints on the characteristics of putative subsurface biospheres.

3.1. The Role of Pressure

With recent advances in high-pressure geomicrobiology, there are grounds for supposing that high pressure may stymie the prospects for life (Picard & Daniel 2013). Hence, it is worth calculating the pressure at the base where T = T_max is attained, and compare it against the known limits for microbes on Earth. The pressure (P_b) at this location is

$$\begin{eqnarray}&&{P}_{b}\approx {\rho }_{c}{g}_{\oplus }{\left(\displaystyle \frac{R}{{R}_{\oplus }}\right)}^{1.7}{H}_{2},\end{eqnarray} \tag{ 13 }$$

where ${\rho }_{c}\approx 2.7\times {10}^{3}\,\mathrm{kg}$ m⁻³ is the roughly constant density of the crust (Whittington et al. 2009), ${g}_{\oplus }$ is the Earth's surface gravity, and H₂ is the lower layer depth of the subsurface habitable region defined in the Appendix; we have also made use of $g\propto M/{R}^{2}$. If we presume that ${H}_{2}\ll R$, we can make use of Equation (8) to obtain

$$\begin{eqnarray}&&{P}_{b}\approx \displaystyle \frac{4\pi {\rho }_{c}{g}_{\oplus }{R}_{\oplus }^{2}F({T}_{\max })}{{\rm{\Gamma }}{Q}_{\oplus }}.\end{eqnarray} \tag{ 14 }$$

An unusual consequence of this equation is the absence of any explicit dependence on R; in other words, to leading order P_b depends only on Γ and T_s; note that the latter has a weak dependence on R as seen from Equation (11).

In principle, one can calculate P_b as a function of T_s, Γ, and R. However, we will focus exclusively on free-floating rocky objects. The reason for doing so is that H₂ in Equation (13) is predicted to be larger, implying that the value of P_b thus obtained constitutes an upper bound for fixed R and Γ. To offer an example, our model yields P_b ≈ 120 MPa for Mars, whereas a free-floating Mars-sized world would have P_b ≈ 460 MPa; in both cases, we have set Γ = 1. Figure 3 depicts P_b as a function of R and Γ. It is apparent that the pressure becomes nearly independent of the radius in accordance with Equation (14). Moreover, for Γ = 1, we find that P_b < 1 GPa holds true for rocky objects with sizes comparable or larger than the Moon.

Zoom In Zoom Out Reset image size

On Earth, the piezophile Thermococcus piezophilus is known to survive at pressures of roughly 125 MPa (Merino et al. 2019), which is comparable to the values in Figure 3 for Γ ≳ 1. Furthermore, laboratory experiments entailing Escherichia coli have found that these bacteria rapidly acquire piezoresistance in the GPa range (Vanlint et al. 2011). Hence, our analysis suggests that the constraints imposed by pressure are insufficient to rule out the prospects for life altogether. One caveat worth noting here is that organisms in deep subsurface environments would need to concomitantly withstand multiple extremes, but our knowledge of polyextremophiles is filled with lacunae (Harrison et al. 2013).

3.2. Biomass in the Deep Biosphere

Ever since the pioneering study by Gold (1992), many studies have sought to quantify the biomass embedded in Earth's subsurface biosphere. The estimates vary by more than an order of magnitude due to the variety of methods deployed for this calculation, as well as the inherent uncertainties surrounding some of the salient parameters (Magnabosco et al. 2018).

A common method for calculating the biomass (M_bio) utilizes the scaling (Gold 1992; Whitman et al. 1998)

$$\begin{eqnarray}&&{M}_{\mathrm{bio}}\propto \displaystyle \frac{{{ \mathcal C }}_{\mathrm{cell}}\,{\delta }_{s}\,{\phi }_{\mathrm{cell}}\,{{ \mathcal V }}_{s}}{{{ \mathcal V }}_{\mathrm{cell}}},\end{eqnarray} \tag{ 15 }$$

where δ_s and ${{ \mathcal V }}_{s}$ represent the porosity and volume of the subsurface habitable region, ϕ_cell denotes the fraction of the pore space occupied by microbes, while ${{ \mathcal C }}_{\mathrm{cell}}$ and ${{ \mathcal V }}_{\mathrm{cell}}$ are the typical organic carbon mass and volume of a typical cell, respectively. It is evident that most of the variables are either geological or biological in nature, and therefore are subject to considerable variability; even on Earth, these parameters are not tightly constrained. The sole exception is ${{ \mathcal V }}_{s}$ because it is possible to loosely estimate it using the results from our model.

Note that ${{ \mathcal V }}_{s}\propto {R}^{2}H$ as it quantifies the area of a spherical shell with thickness H. Although H is a complex multivariable function, we can draw upon the preceding results. For T_s < T_min—which is true for most rocky objects of interest herein—it was shown in Figure 1 that H becomes effectively independent of T_s. This leaves us with the parameters R and Γ. From Equations (8) and (10), which are fairly accurate for $R\gtrsim 0.1\,{R}_{\oplus }$, we see that the same scaling applies, namely, $H\propto {R}^{2}/Q$.⁴ Hence, from these scalings, we find that ${{ \mathcal V }}_{s}\propto {R}^{4}/Q\propto {R}^{0.3}/{\rm{\Gamma }}$ with the last relationship following from Equation (4).

Thus, with all other factors in Equation (15) held fixed apart from ${{ \mathcal V }}_{s}$, the biomass is estimated as

$$\begin{eqnarray}&&{M}_{\mathrm{bio}}\sim 20\,\mathrm{Pg}\,{\rm{C}}\left(\displaystyle \frac{1}{{\rm{\Gamma }}}\right){\left(\displaystyle \frac{R}{{R}_{\oplus }}\right)}^{0.3},\end{eqnarray} \tag{ 16 }$$

where the normalization was adopted from the mean value of Magnabosco et al. (2018, pg. 712). A striking aspect of Equation (16) is the relatively weak dependence on R. In principle, a Super-Earth with R ≈ 1.5 R_⊕and Γ ≈ 0.2 might support a subsurface biosphere whose biomass is ∼20% that of Earth's total biomass, which equals 550 Pg C (Bar-On et al. 2018). In actuality, however, we caution that M_bio is constrained by energy, nutrient, and miscellaneous limitations, owing to which M_bio constitutes an idealized loose upper bound.

Neither energy nor nutrient limitation are easy to model, and require a case-by-case treatment. Instead, we will focus on a single energetic pathway. Among the various energy sources, radiolysis has garnered widespread attention (Lin et al. 2005; Sherwood Lollar et al. 2007; Atri 2016; Dzaugis et al. 2018), especially after the discovery of Desulforudis audaxviator, a sulfate-reducing bacterium living 2.8 km beneath Earth's surface that derives its energy from this source (Chivian et al. 2008). The maximum subsurface biomass (M_max) that might be supported by methanogenesis with the reactant (H₂) derived from radiolysis is

$$\begin{eqnarray}&&{M}_{\max }\approx \displaystyle \frac{{{ \mathcal C }}_{\mathrm{cell}}\,{\rm{\Delta }}G\,\bar{{ \mathcal P }}\,{{ \mathcal V }}_{s}}{{P}_{\mathrm{bio}}},\end{eqnarray} \tag{ 17 }$$

where $\bar{{ \mathcal P }}$ is the average volumetric production rate of H₂ via radiolysis, ΔG is the Gibbs free energy of methanogenesis, and P_bio is the minimal power required for the viability of a single microbe. As most of these parameters are unknown for other objects, we hold them fixed at the average values on Mars and/or Earth. In the Noachian era, $\bar{{ \mathcal P }}$ was probably ≲10 times higher than today due to the higher abundance of radionuclides (Tarnas et al. 2018), while the other parameters were ostensibly similar. Thus, by this reckoning, the total H₂ that might arise from radiolysis on modern Mars is roughly an order of magnitude lower than the Noachian mean value of ∼3 × 10¹⁰ moles yr⁻¹ (Tarnas et al. 2018). It was estimated by Sholes et al. (2019) that a total H₂ flux of ∼1.3 × 10⁸ moles yr⁻¹ could support ∼10²⁷ cells. By employing these predictions along with ${{ \mathcal C }}_{\mathrm{cell}}\sim 2.1\times {10}^{-17}$ kg (Magnabosco et al. 2018, pg. 707), $\bar{{ \mathcal P }}\propto {\rm{\Gamma }}$ (because the production rate is proportional to radionuclide concentration) and the scaling ${{ \mathcal V }}_{s}\propto {R}^{0.3}/{\rm{\Gamma }}$ derived earlier, we arrive at

$$\begin{eqnarray}&&{M}_{\max }\sim 0.6\,\mathrm{Pg}\,{\rm{C}}{\left(\displaystyle \frac{R}{{R}_{\oplus }}\right)}^{0.3}.\end{eqnarray} \tag{ 18 }$$

Upon comparing Equations (18) and (16), we find that M_max is ∼3% of Earth's subsurface biosphere for $R\approx {R}_{\oplus }$. Moreover, at this radius, M_max is merely ∼10⁻³ that of Earth's global biomass (Bar-On et al. 2018). The weak scaling in Equation (18) implies that M_max is only ≲2 times smaller on objects such as the Moon ($R=0.27\,{R}_{\oplus }$) and Mars ($R=0.53\,{R}_{\oplus }$). Hence, in principle, both these objects might have been capable, at some point perhaps, of hosting deep biospheres that were approximately three orders of magnitude smaller than Earth's global biomass.

In closing, we reiterate that only energetic and thermal (which dictates ${{ \mathcal V }}_{s}$) constraints were considered. Aside from limitations imposed by nutrients and reactants, liquid water is obviously another vital requirement. However, this issue is hard to quantify because it depends, inter alia, on the pathways of object formation, volatile delivery by impactors, crust composition, and potentially the deep water cycle (Peslier et al. 2017). It is estimated that Earth's crust comprises ∼4 × 10²³ mL of water (Peslier et al. 2017, their Table 1). If roughly ∼10% of water is present in the subsurface biosphere region and it hosts ∼10³⁰ cells, we find that the average cell density is ∼2.5 × 10⁷ cells/mL; this estimate is compatible with sophisticated models (Magnabosco et al. 2018). Hence, provided that other rocky objects are not significantly devoid of volatiles, it seems unlikely prima facie that subsurface biospheres would be limited by water availability.

In this work, we examined the prospects for life in deep terrestrial environments. While this subject has been investigated for specific rocky objects (e.g., Earth), we constructed simple models to assess the salient features of such putative biospheres.

We began by studying the extent of the habitable region (H) associated with deep biospheres as a function of the surface temperature, relative radionuclide abundance (Γ), and size (R) of the world. We showed that H may span orders of magnitude for different combinations of these parameters. For rocky objects that are Moon-sized and larger, we determined that the scaling $H\propto {R}^{-1.7}{{\rm{\Gamma }}}^{-1}$ is valid. We compared the results from our model with prior studies (involving Mars) and found that they exhibit good agreement. Then we investigated two other constraints on the habitability of these settings. First, we studied the pressure at the base of the biosphere and showed that it is nearly independent of R. We determined that pressure is unlikely to completely rule out the prospects for putative lifeforms. Lastly, we presented simple estimates for the total biomass that could be supported in deep terrestrial habitats. We determined that this biomass may depend weakly on R, and that it might be a few percent of Earth's deep biosphere and three orders of magnitude smaller than Earth's global biomass.

In closing, two factors merit highlighting. First, the number of rocky planets in the HZs of host stars is ≲0.3 per star (Zink & Hansen 2019). Even if not all of these objects are capable of sustaining liquid water on the surface over long timescales (e.g., Mars), they may still have the capacity to support deep biospheres. Furthermore, gravitational microlensing studies (Strigari et al. 2012; Bhatiani et al. 2019) and extrapolation of the size distribution function of interstellar objects (Rice & Laughlin 2019; Siraj & Loeb 2019) both indicate that the number of free-floating objects with sizes larger than the Moon is $\gtrsim 10$-100 per host star.⁵ Even if only a small fraction of them were formed inside the snow line (prior to ejection) and are predominantly rocky, they may nevertheless outnumber planets in the HZs of stars.

Second, because deep biospheres are situated underneath the surface, detecting unambiguous signatures of biological activity is not readily feasible via remote sensing techniques. The most obvious solution is to carry out in situ studies of rocky objects in our backyard such as the Moon and Mars. The Moon was habitable shortly after its formation (Schulze-Makuch & Crawford 2018), and it is not altogether inconceivable that some traces and markers of life might survive to this day. NASA's Artemis program can pave the way for lunar subsurface exploration via establishing a sustainable base on the Moon purportedly by 2024.⁶

While methods like transient electromagnetic sounding can reveal the presence of subsurface water, advances in drilling technologies are necessary to develop autonomous machines capable of reaching km-scale depths (Carrier et al. 2020). Alternatively, one may carry out excavations at the sites of large craters on Mars to extract and analyze materials derived from the subsurface (Cockell & Barlow 2002). A more ambitious (and technically feasible) enterprise entails combing the Oort cloud for the nearest Moon-sized extrasolar objects and sending fast spacecraft with appropriate payloads to scrutinize them (Hein et al. 2020).

This work was supported by the Breakthrough Prize Foundation, Harvard University's Faculty of Arts and Sciences, and the Institute for Theory and Computation (ITC) at Harvard University.

As noted in Section 2, there are two cases that we need to consider in order to calculate the spatial extent of the subsurface habitable region (H). We examine them in detail below.

In the first case, T_s > T_min. In this context, the top layer of the habitable region starts at the surface, i.e., at r = R. Let us suppose that the bottom layer of the habitable region—where T = T_max is attained—occurs at a depth H beneath the surface (namely r = R − H). In this particular scenario, therefore, the spatial extent of the habitable subsurface region is H and is found by setting T = T_max and r = R − H in Equation (6) and using the definition of F(T) introduced just below (6). After doing so, we end up with

$$\begin{eqnarray}&&\displaystyle \frac{Q\left[{R}^{2}-{(R-H)}^{2}\right]}{8\pi {R}^{3}}=F({T}_{\max }),\end{eqnarray} \tag{ A1 }$$

which can be further simplified to yield

$$\begin{eqnarray}&&\displaystyle \frac{Q\left(2{RH}-{H}^{2}\right)}{8\pi {R}^{3}}=F({T}_{\max }),\end{eqnarray} \tag{ A2 }$$

which is a quadratic equation in H. By solving this equation for H and eliminating the non-physical root (which exhibits H > R), we end up with Equation (7).

Now, let us tackle the second setting wherein T_s < T_min. In this case, the top layer of the habitable subsurface region will not be at the surface, but will begin at a depth H₁ (namely r = R − H₁) below it. One can solve for H₁ as follows: at the location r = R − H₁, the temperature is equal to T_min. By substituting this value of r and T in Equation (6), and repeating the same procedure in the above paragraph, we find

$$\begin{eqnarray}&&\displaystyle \frac{{H}_{1}}{R}=1-\sqrt{1-\displaystyle \frac{8\pi R}{Q}F({T}_{\min })}.\end{eqnarray} \tag{ A3 }$$

Next, the bottom layer of the habitable surface region will occur when the temperature is equal to T_max. If we denote the depth below the surface where this occurs by H₂ (corresponding to $r=R-{H}_{2}$), we can solve for the latter quantity by demanding that T = T_max at r = R − H₂. After making use of Equation (6) and carrying out the resultant algebra, we end up with

$$\begin{eqnarray}&&\displaystyle \frac{{H}_{2}}{R}=1-\sqrt{1-\displaystyle \frac{8\pi R}{Q}F({T}_{\max })}.\end{eqnarray} \tag{ A4 }$$

In this scenario, recall that the habitable subsurface region begins at H₁ and is terminated at H₂. Therefore, the spatial extent of the subsurface habitable region is H = H₂ − H₁, which leads us to Equation (9).

For free-floating worlds, the second case is indeed valid. Here, let us recall that

$$\begin{eqnarray}&&F({T}_{\min })={ \mathcal A }\left({T}_{\min }-{T}_{s}\right)+{ \mathcal B }\mathrm{ln}\left(\displaystyle \frac{{T}_{\min }}{{T}_{s}}\right),\end{eqnarray} \tag{ A5 }$$

$$\begin{eqnarray}&&F({T}_{\max })={ \mathcal A }\left({T}_{\max }-{T}_{s}\right)+{ \mathcal B }\mathrm{ln}\left(\displaystyle \frac{{T}_{\max }}{{T}_{s}}\right),\end{eqnarray} \tag{ A6 }$$

where T_s in both of the above equations is set by Equation (11). One can then substitute these two expressions into Equation (9) to solve for H as a function of R and Γ. Due to the cumbersome nature of the expression, it is not explicitly provided here.