© ESO 2018

1 Introduction

The recent discovery of seven roughly Earth-sized planets orbiting the low-mass star TRAPPIST-1 (Gillon et al. 2017) has vaulted this system to the forefront of exoplanetary characterization. The planets orbit the star with semi-major axes <0.1 AU, and orbital periods of a few Earth days. Given their proximity to the star, and the star’s low mass and luminosity, the surface of each planet has a moderate temperature, ranging from ~160 to 400 K (Gillon et al. 2017). Given the planets’ mean densities, these low temperatures suggest that some might have solid surfaces composed of H₂O ice and/or rock. The planets also have non-zero orbital eccentricities (Gillon et al. 2017; Wang et al. 2017), so that the tidal forces on the bodies vary with time, resulting in heating of their interiors by tidal dissipation (Luger et al. 2017). Tidal heating is an important energy source in the satellites of the outer planets in our solar system (e.g., McEwen et al. 2004; Greeley et al. 2004; Spencer et al. 2009) and could significantly enhance the habitable volume in the galaxy by warming exoplanets and their satellites (exomoons) (Dobos & Turner 2015).

Here, tidal dissipation (or tidal heating) refers to the process of the dissipation of orbital energy in the interior of a solid body (a moon or a planet). For a single secondary object, that is, a lone moon orbiting a planet, or a lone planet orbiting a star, tidal dissipation rapidly decreases the orbital eccentricity of the secondary, and then ceases when e = 0 (Murray & Dermott 1999).

If multiple objects orbit the primary, and the orbiting objects occupy a mean motion resonance, periodic gravitational perturbations will help maintain non-zero eccentricities. This is the case, for example, for the satellites of Jupiter and Saturn. Jupiter’s innermost moon, Io, has a tidal heat flux about ~1–2 W m⁻² (Spencer et al. 2000; Veeder et al. 2004). The heat flux is so large that the moon has become volcanically active, erupting silicate magmas and sulfur-rich compounds onto its surface (McEwen et al. 2004). Tidal heating can also be an important heat source in bodies with a significant surface layer of H₂ O (whether it be solid ice or liquid water). Saturn’s icy moon Enceladus, which has a mean radius of 252 km (Porco et al. 2006) and is thus only asteroid-sized, has plumes of water-rich material erupting from tidally induced cracks at its south pole (Porco et al. 2006), and a localized heat flow between 0.1 and 0.25 W m⁻², roughly three times that of the Earth (Spencer et al. 2009). Tidal deformation and the geodynamical consequences of tidal heating have also created the bizarre geology on the surfaces of Jupiter’s moons Europa (Greeley et al. 2004) and Ganymede (Pappalardo et al. 2004).

The TRAPPIST-1 planets occupy a mean motion resonance chain with orbital period ratios slightly larger than 24:15:9:6:4:3:2 (Unterborn et al. 2018; Luger et al. 2017). The resonances maintain the orbital eccentricities of the planets resulting in long-term tidal heating (Luger et al. 2017). Luger et al. (2017) find that assuming Earth-like tidal dissipation rates, TRAPPIST-1b may have a tidal heat flux larger than Io (~ 3 W m^-2) with peaks larger than 10 W m^-2 due to the periodic change of its orbital eccentricity. The other six planets in the system would experience weaker tidal heating, with heat fluxes ranging between values similar to the Earth and Io for planets c–e, and less than Earth for planets f–h.

Here, we constructed simple compositional and thermal models of the interiors of each of the TRAPPIST-1 planets. Assuming that the planets are composed of uniform-density (incompressible) rock, iron, and H₂ O, we explored the range of possible compositional and interior models permitted by the error bars in mass and radius for each planet. Because the error bars on mass and radius are so large, more sophisticated models including, for example, compression and thermalevolution, are unlikely to give robust and meaningful results. Our calculations allowed us to determine the precision with which the mass of each planet must be measured to determine whether a given planet harbors a significant reservoir of H₂ O. We computed the amount of tidal dissipation in each object based on a uniform viscosity and rigidity, which depend on the volume fractions of ice and rock in each planet. We also determined the expected temperature in the rock mantles of planets b, c, d, e, and f by setting the tidal heat flux equal to the amount of heat that could be removed by solid-state convection. Some of the planets have globally averaged heat fluxes greater than the heat flux that can trigger a runaway greenhouse effect, which causes irreversible water evaporation and loss (Kasting et al. 1993).

2 Interior structures

2.1 Prior work

Table 1 summarizes the values of mass, radius, mean density, effective temperature, and orbital period for each of the seven TRAPPIST-1 planets. The vast majority of the values originate from Wang et al. (2017). The masses range from almost 0.1 to 1.63 M_⊕, where M_⊕ = 5.98 × 10²⁴ kg is the mass of the Earth. The radii are generally close to the radius of the Earth (R_⊕ = 6371 km).

Several studies have constrained interior states for the TRAPPIST-1 planets using PREM (Preliminary Earth Reference Model), a density and seismic velocity profile for the Earth (Dziewonski & Anderson 1981). Zeng et al. (2016) use an equation of state for silicate material based on PREM to constrain iron contents and core sizes in rocky exoplanets as a function of their size and temperature. These results have been used by Gillon et al. (2017) and Wang et al. (2017) to constrain the rock/metal fractions in the TRAPPIST-1 planets. Wang et al. (2017) combine the original discovery data for the TRAPPIST-1 system with 73.6 days of K2 observations to obtain updated masses and radii for the planets, which serve as the primary source of mass and radius information for our study. Wang et al. (2017) also use equations of state for rock and metal based on PREM to construct interior models of the planets, demonstrating that several of the bodies have densities low enough to contain a significant amount of water in their interiors.

Another method is to use equations of state for silicate minerals (e.g., MgSiO₃) to calculate density and pressure as a function of depth, constrained by the mass and radius of the planet (Unterborn et al. 2018). This detailed analysis has been performed for TRAPPIST-1b and 1c. In these models, phase transitions from one polymorph of MgSiO₃ to another are considered, resulting in a mantle which contains layers of enstatite, pyroxene, and garnet, before a final transition to Bridgmanite at a depth of ~ 1000 km. Each of these phaseshave much lower densities than compressed Bridgmanite, leading to a low-density upper mantle, and permitting a relatively thin H₂O layer.

2.2 Methods

Here, we have calculated the simplest possible interior structures for each of the TRAPPIST-1 planets assuming that the rock, metal, and ice layers in their interiors have uniform densities. We calculated the suite of structures that fit the mean mass and radius, and explore the range of possible interior structures based on the uncertainties of planetary masses and radii. Given the large error bars on the masses and radii of the planets, some of which permit planetary mass of zero, simple uniform-density models provide interior structures that are just as valid as those takinginto account the pressure and temperature changes with depth inside the bodies (see e.g., Unterborn et al. 2018).

Table 1 shows the mean densities for the planets that we calculated based on their maximum and minimum permitted masses and radii; the large error bars on these quantities yield large error bars on the mean densities, as well. Despite the uncertainties, the densities of the TRAPPIST-1 planets show quite a bit of variation, ranging from values close to the density of water ice under compression, to densities hinting at rocky, and possibly iron-rich compositions.

Our choices of densities for rock, iron, and ice include the effects of compression and phase changes that will occur at depth inside Earth-sized planets. For iron, we assumed a constant density of ρ_Fe = 12 000 kg m⁻³, consistent with the density of Earth’s inner core (Dziewonski & Anderson 1981; Unterborn et al. 2018). For rock in planets b, c, d, and g we assumed a constant density of ρ_r = 5000 kg m⁻³, appropriate for compressed Bridgmanite (MgSiO₃) (McCammon et al. 2016; Unterborn et al. 2018). For H₂ O at low pressures, we assumed a density of ρ_I = 1000 kg m⁻³, which accounts for the possibility of both a solid ice I shell and a liquid water ocean. For H₂ O at high pressures, above 209 MPa, we assumed a density of ρ_hpp = 1300 kg m⁻³, which accounts for the possible presence of high-pressure ice polymorphs (hpp) II through VII (Hobbs 1974), each of which could be present in distinct layers within the planet. When high-pressure ice phases are present, it is likely that ices VI and VII will dominate because they are stable over a broad range of temperatures at the ~ tens of GPa central pressures appropriate for TRAPPIST-1 bodies.

For planets e, f, and h, the minimum permitted planetary densities based on the masses and radii from Wang et al. (2017) are smaller than the density of high-pressure H₂ O polymorphs. In the context of our assumptions, a planet with $\overline{\rho }={\rho }_{\text{hpp}}$ would correspond to a planet composed of pure H₂O, which we consider unlikely based on planet and satellite formation models (e.g., Canup & Ward 2002; Morbidelli et al. 2012; Johansen et al. 2015). In Table 1, we reported a minimum mass corresponding to the mass of a planet with the maximum allowed radius, and $\overline{\rho }\approx {\rho }_{\text{hpp}}$.

When determining the interior structures in the low-density limit for planets e and f, and for all of the interior structures for planet h, we adopted a density for rock ρ_r = 3300 kg m⁻³, the density of Prinn–Fegley rock. Prinn–Fegley rock has been used as a representative composition for the rocky component of the outer planet satellites (Mueller & McKinnon 1988; Barr & Canup 2008). This could represent a more hydrated, but compressed rock core.

With the exception of TRAPPIST-1c, all of the seven planets have densities low enough to indicate the presence of H₂ O (Wang et al. 2017; Unterborn et al. 2018). For these bodies, we began by estimating the maximum depth of the ice I/liquid water layer: z_I = P_I∕(ρ_Ig), where $g=GM/R_{\text{m}}^2$ is the surface gravity, and P_I = 209 MPa is the maximum depth at which ice I is possible (Hobbs 1974). Subtracting the mass of ice I, we determined the mean density of the remaining planetary material, which allowed us to constrain the relative sizes of the mantles of high pressure ice polymorphs and rock, and the radius of the iron core. Interior models applied by Wang et al. (2017) suggest that TRAPPIST-1b may contain approximately 25% water by mass (with a large uncertainty between about 0% and 60%). Using this constraint of 25% mass fraction, we were able to obtain a unique interior structure for TRAPPIST-1b that satisfied the mass and radius constraints.

3 Thermal model

Here, we have computed the heat production by tides and heat lost by conduction and convection to constrain tidal heating rates and rock mantle temperatures in each of TRAPPIST-1 exoplanets. The most common method of estimating tidal heating is to assume that the interior of the secondary object (in this case, the planet) has a Maxwell viscoelastic rheology (Love 1906; Peale & Cassen 1978), characterized by a shear modulus, μ, and viscosity, η. More complex rheologies (for example, Andrade or a standard linear solid) can also be used, and yield slightly different behaviors (e.g., Efroimsky 2012). However, given the large uncertainties in the planets’ compositions, we chose the Maxwell model because it is a simple model with simple behavior and relatively few parameters.

In a Maxwell viscoelastic material, heating is maximized in material whose Maxwell time, τ = η∕μ is comparable to the orbital period of the secondary, P (Love 1906). For ice I, the low-density phase of H₂O that floats on liquid water, the Maxwell time is comparable to the orbital period of Europa, 85 h, if the ice is close to its melting point (Barr & Showman 2009).

In a tidally heated body, the interior temperature depends upon the balance and feedback between internal heat generation by, for example, radiogenic and tidal sources, and heat transported by conduction and possibly solid-state convection (Schubert et al. 2001; Tackley 2001; Hussmann et al. 2002; Moore 2006). A body whose interior, or portion of its interior, has a viscosity and rigidity close to the values for which P ~ τ, will undergo intense tidal heating. A body with a slightly higher viscosity, that is, one that is cooler, can be warmed by tidal heating until P ~ τ. The body can experience inner melting, which in silicate bodies leads to volcanism (Peale 2003), and in ice/rock bodies, leads to the maintenance of liquid water oceans (Ojakangas & Stevenson 1989). Thus, the tidal heating rate in a solid body and its interior structure (namely, composition, temperature, and phase as a function of depth), depend on one another, and the body may reach an equilibrium structure where heat generation is equal to heat loss.

Figure 1 shows a schematic illustration of the methods we used to calculate the internal temperature inside each of the TRAPPIST-1 planets based on their different interior structures. Given values for the volume fractions of each of the planets’ constituents (ϕ), we constructed a volume-averaged viscosity (η) and rigidity (μ). Based on these average values, we calculated the response of the planet to the tidal forcing, expressed by the imaginary part of the k₂ Love number, Im(k₂). This value was used to calculate the surface heat flux from tidal dissipation F_tidal. The heat generated from tidal dissipation is likely removed from the planets’ deep interiors via solid-state convection, the dominant heat transfer mechanism in the terrestrial planets and ice/rock satellites in our solar system. To determine the temperature in the planets’ deep interiors, we balanced F_tidal and the convective heat flux (F_conv), and found the value of temperature (T_eq) for which F_conv = F_tidal. If T_eq is above thesolidus temperature for rock, the planet may possess a magma ocean in its interior. The globally averaged energy flux on the planet is the sum of the energy flux from the mean stellar irradiation (F_⋆) and tidal heating, $${\overline{F}}_{\text{glob}}=F_{\star }+F_{\text{tidal}},$$(1)

where F_⋆ is taken at the top of the planet’s atmosphere. The stellar radiation can be expressed as $$F_{\star }=\frac{L_{\star }\left(1-A_{\text{B}}\right)}{16\pi a^2\sqrt{1-e^2}},$$(2)

where L_⋆ is the stellar luminosity, A_B is the Bond albedo of the planet, a and e are the semi-major axis and the eccentricity of the planet’s orbit, respectively. Greenhouse and cloud feedbacks were neglected in the calculation of F_⋆. The albedo of each of the planets was set to 0.3, an Earth-like value.

Fig. 1 Schematic representation of our tidal and thermal model. Left: internal structure for the planet, composed of water ice I (white), high-pressure ice polymorphs (HPP, gray), rock (dark gray), and iron (black). The volume fractions of each constituent material were calculated based on the mean density of each planet (see Sect. 2). Middle: to calculate the response of each planet to tidal forcing, we constructed a volume-weighted average viscosity and rigidity for the planet based on the ϕ values and the rheology of each constituent material. A value of the imaginary part of the k2 Love number, and the tidal heat flux Ftidal were calculated. Right: the temperature in the planet’s interior (Teq) was governed by a balance between heat generation by tides and transport by solid-state convection; the efficiency of both processes depends sensitively on the viscosity and rigidity of the planet. The net global heat flux ${\overline{F}}_{\text{glob}}$ was calculated from the tidal heat flux and the incident sunlight from the star, F⋆.

3.1 Heat generation

Here, we have focused on tidal heating as the dominant source of energy for heating the interiors of the TRAPPIST-1 planets. This is because the magnitude of tidal heating expected in rock- or ice-dominated bodies at the orbital periods associated with the TRAPPIST-1 planets, a few to ten watts per meter squared, is several orders of magnitude higher than the heat flows associated with radiogenic heating, which is typically of the order of tens of milliwatts per meter squared (Henning et al. 2009; Henning & Hurford 2014). Residual heat of accretion is also a possible heat source in any solid planet. The energy per unit mass released during accretion, E_acc ~ (3∕5)GM∕R is sufficient to melt planets as large as the TRAPPIST-1 bodies, and is still a significant contribution to the Earth’s heat flux. Even if the TRAPPIST-1 bodies still retained 20% of their accretional heat, the resulting heat flux would still be much smaller than the tidal heat flux.

The tidal heat flux was calculated using a viscoelastic model for a homogeneous body described by Dobos & Turner (2015), which was originally developed by Henning et al. (2009). Although the planet was assumed to be homogeneous, we mimicked the effect of multiple materials by determining the effective viscosity and rigidity for each layers and weighting them by their volume fractions. To approximate the material properties of a planet composed of several different materials, we calculated the volume fraction of each material contained in the planet: for a given material, i, the volume fraction ϕ = V_i∕V _tot, where $V_{\text{tot}}=(4/3)\pi R_{\text{pl}}^3$ is the total volume of the planet. We calculated values of the volume fractions of ice I (ϕ_I), high-pressure ice polymorphs (ϕ_hpp), rock (ϕ_r), and iron (ϕ_Fe) contained ineach planet, using the mass and radius of the planet as constraints, assuming uniform densities for each of the materials. Asingle uniform viscosity and rigidity for the planet was approximated by $$\eta \approx {\varphi }_{\text{I}}{\eta }_{\text{I}}+{\varphi }_{\text{hpp}}{\eta }_{\text{hpp}}+{\varphi }_{\text{r}}{\eta }_{\text{r}}+{\varphi }_{\text{Fe}}{\eta }_{\text{Fe}},$$(3)

where η is the viscosity of the material. We used a similar relationship to construct a single value of the shear modulus (μ) that approximates the behavior of the entire planet: $$\mu \approx {\varphi }_{\text{I}}{\mu }_{\text{I}}+{\varphi }_{\text{hpp}}{\mu }_{\text{hpp}}+{\varphi }_{\text{r}}{\mu }_{\text{r}}+{\varphi }_{\text{Fe}}{\mu }_{\text{Fe}}.$$(4)

The total amount of energy produced by tidal dissipation in a synchronously rotating body (Segatz et al. 1988), $${\dot{E}}_{\text{tidal}}=-\frac{21}{2}\text{Im}(k_2)\frac{R^5{\omega }^5{e}^2}{G},$$(5)

where ω = 2π∕P is the orbital frequency of the planet, P is the planet’s orbital period, G is the gravitational constant, e is the planet’s eccentricity, and R is its radius. The quantity Im(k₂) is the imaginary part of the k₂ Love number, which describes how the planet’s gravitational potential is disturbed by its tidal distortion. The value of Im (k₂) depends on the structure and rheology in the body (Segatz et al. 1988).

There are a variety of different ways of calculating Im(k₂) for a solid planet. In some studies, Im(k₂) is expressed as k₂∕Q, where Q is a constant parameter describing the fraction of orbital energy dissipated per tidal cycle (see e.g., Murray & Dermott 1999 for discussion, and Papaloizou et al. (2017) for example application of this approach to the TRAPPIST-1 system). The value of Q for solid planetary bodies is not well-constrained, and is generally thought to be between Q ~ 10 and 200 (Goldreich & Soter 1966; Murray & Dermott 1999).

A constant-Q approach, however, does not take into account the feedback between the thermal state of a planet and its response to tidal forces. A warm planet will deform more in response to tidal forces than a cold, rigid planet. Thus, it is more accurate to calculate Im (k₂) based on the viscosity and rigidity of the planet’s interior. A common approach is to assume that the planet is composed of a single material (for example, rock). Given the density, shear modulus, and viscosity of the planetary material, and an assumption about the rheology (for example Maxwell, Andrade, etc., Efroimsky (2012)), it is possible to derive an analytic expression for Im (k₂) (e.g., Harrison 1963; Henning et al. 2009; Dobos & Turner 2015, see Eq. (6) here).

A more sophisticated approach is to assume that the planet is composed of, perhaps two or three materials, each of which has a different density, shear modulus, and viscosity. In this case, analytic, although quite cumbersome expressions can still be obtained for Im (k₂) (Harrison 1963; Remus et al. 2012, 2015; Beuthe 2013). Many of these studies rely upon simplifying assumptions to obtain closed-form analytic solutions (e.g., Remus et al. 2012, 2015); in most cases, these simplifying assumptions may not apply to the TRAPPIST-1 planets. Another approach is to calculate Im (k₂) numerically, by representing the planet by ten or more layers, each of which has a density, shear modulus, and viscosity dictated by the composition of the layer, but also thermal state of the material (e.g., Sabadini et al. 1982; Moore 2006; Roberts & Nimmo 2008).

Given the large uncertainties in the masses and radii of the TRAPPIST-1 planets, calculation of interior thermal states using more sophisticated models (for example, one in which the heat generation and transport in each different compositional layer inside the planet is treated separately) would be premature. Therefore, we used a Maxwell viscoelastic rheology to calculate the value of Im (k₂): $$-\text{Im}(k_2)=\frac{57\eta \omega }{4\overline{\rho }gR\left[1+{\left(1+\frac{19\mu }{2\overline{\rho }gR}\right)}^2\frac{{\eta }^2{\omega }^2}{{\mu }^2}\right]},$$(6)

where g = GM∕R² is the planet’s surface gravity, η is the viscosity (Eq. (3)) and μ (Eq. (4)) is the shear modulus of the planet (Henning et al. 2009).

The globally averaged tidal heat flux was calculated from $$F_{\text{tidal}}=\frac{{\dot{E}}_{\text{tidal}}}{4\pi R^2}.$$(7)

It is worth noting that this tidal heating model is valid only for small orbital eccentricities, that is, for e ≲ 0.1. At higher eccentricities, additional higher-order terms must be included in the tidal potential, which can increase the amount of heating (Mignard 1980). Fortunately, this constraint is fulfilled for each of the seven planets, according to the data of Gillon et al. (2017) and Wang et al. (2017). In our model, we did not consider gravitational interactions of the planets with each other, only the star-planet (two-body) tidal interaction was considered in each case.

3.2 Heat loss

The tidal heat generated inside the planet was assumed to be transported to the surface by solid-state convection, then conducted across a surface boundary layer where it is radiated into space. The heat flux across a convecting mantle, $$F_{\text{conv}}=\frac{k_{\text{therm}}\Delta T}{D}\text{Nu},$$(8)

where k_therm is the thermal conductivity, Δ T = T_mantle − T_surf is the difference in temperature between the hot, convecting mantle (T_mantle) and T_surf is the temperature at the surface, and D is the thickness of the convecting layer. The parameter Nu, the Nusselt number, describes the relative efficiency of convective heat transport versus conductive heat transport. In an internally heated planet with a Newtonianrheology, Nu is related to the thermal, physical, and rheological properties of the mantle (Solomatov & Moresi 2000), $$\text{Nu}=0.53{\theta }^{-4/3}{\text{Ra}}^{1/3},$$(9)

where $\theta =(Q^{*}\Delta T/R_{\text{G}}{T}_{\text{mantle}}^2)$, Q^* is the activation energy in the flow law for the convecting material and R_G = 8.314 J⁻¹ (mol K)⁻¹ is the universal gas constant. The parameter Ra is the Rayleigh number, which describes the vigor of convection, $$\text{Ra}=\frac{\rho g\alpha \Delta T{D}^3}{\kappa \eta (T_{\text{mantle}})},$$(10)

where the coefficient of thermal expansion α, thermal diffusivity κ_therm = k_therm∕(ρC_p), C_p, and η(T_mantle) is the viscosity evaluated at the mantle temperature. Using Eq. (10) to evaluate the Rayleigh number in Eq. (9), we obtained an expression for the convective heat flux which did not depend on the thickness of the convecting layer, or the temperature at the surface of the planet (Solomatov & Moresi 2000; Barr 2008), $$F_{\text{conv}}=0.53{(\frac{Q^{*}}{R_{\text{G}}{T}_{\text{mantle}}^2})}^{-4/3}{(\frac{\rho g\alpha k_{\text{therm}}^3}{{\kappa }_{\text{therm}}\eta (T_{\text{mantle}})})}^{1/3}.$$(11)

We assumed that the majority of tidal heat will be produced in, and transported within, the rock mantle of the planet. Thus, the parameters in Eq. (11) were evaluated using the properties of rock: Q^* = 333 kJ mol⁻¹ (see Sect. 3.4.1), k_therm = 3.2 W m⁻¹ K⁻¹, α = 3 × 10⁻⁵ K⁻¹, and C_p = 1200 J kg⁻¹ (Solomatov & Moresi 2000).

To determine the equilibrium (or equilibria) between F_conv and F_tidal, we evaluated Eqs. (7) and (11) for a variety of values of T_mantle, the temperature in the rock mantle of each planet. Two equilibria are possible between the tidal heat flux and convective cooling flux which correspond to two equilibrium temperatures of T_mantle, which we referred to as T_eq (Moore 2006; Dobos & Turner 2015). The left panel in Fig. 2 illustrates schematically how the tidal and convective heat fluxes vary as a function of temperature. One equilibrium can exist at a temperature well below the solidus, corresponding to a balance between tidal heat generation and transport in a purely solid planetary mantle (point (A) in the left panel, and also see the middle panel in Fig. 2). In this case, the amount of tidal heat generated is sufficiently low to be removed by convection, even if the convective flow is relatively sluggish. This equilibrium is unstable – if the tidal heat generation increases, convection cannot remove the additional tidal heat (Moore 2003; Dobos & Turner 2015). Depending on the tidal forcing frequency and the rheology of the planet, a second equilibrium can be achieved in which the planet is partially molten (point (B) in the left panel, and the right panel in Fig. 2). Tidal heat generation decreases as a function of temperatureand melt fraction above the solidus because partially molten rock is less dissipative than warm, solid rock. However, the convective heat flux increases sharply as the presence of melt decreases the mantle viscosity. This high-temperature equilibrium is stable (Moore 2003; Dobos & Turner 2015) – as the tidal heating rate increases, the viscosity of the mantle decreases, permitting more efficient convective heat flow and resulting in cooling of the mantle.

3.3 Runaway greenhouse limit

To characterize the habitability of the planets, we determined whether the planet might enter a runaway greenhouse state. When the globally averaged heat flux on a planet exceeds a critical value, then surface water slowly evaporates. The water molecules inthe upper atmosphere disintegrate to hydrogen and oxygen due to photodissociation by ultraviolet photons, and then hydrogen molecules can easily escape to space. Hence, entering the runaway greenhouse phase could imply that water molecules would not be able to re-form later, even if the temperature decreases (Kasting et al. 1993).

To determine the runaway greenhouse limit, we used the formulation of Pierrehumbert (2010, Chap. 4). Ifthe globally averaged energy flux on the planet from stellar irradiation and tidal heat (${\overline{F}}_{\text{glob}}$) exceeded the minimum energy flux on top of a water-rich atmosphere then it triggers the runaway greenhouse effect: $$F_{\text{RG}}=o\frac{\sigma {(l/R_{\text{water}})}^4}{\ln {\left(p^{\star }/\sqrt{2p_{0}g/{\kappa }_0}\right)}^4},$$(12)

where o = 0.7344, σ is the Stefan–Boltzmann constant, l = 2.425 × 10⁶ J kg⁻¹ is the latent heat, R_water = 461.5 J (kg K)⁻¹ is the gas constant for water vapor, κ₀ = 0.055 is the gray absorption coefficient, p₀ = 10 000 Pa is the reference pressure for absorption, $p^{\star }=p_{\text{ref}}\exp \left(\frac{l}{R{T}_{\text{ref}}}\right)$, p_ref = 610.616 Pa and T_ref = 273.13 K.

Fig. 2 Left: schematic representation of the heat flux (in log units) from tidal heating (black line) and convective heat transport (gray) as a function of temperature. Two equilibria are possible: point (A) represents a balance between tidal heating and convection in a purely solid mantle (Teq is less than the solidus temperature Ts). This stateis depicted schematically in the middle panel. Point (B) represents an equilibrium for which Teq > Ts, implying a partially molten rock mantle exists within the planet. This state is shown schematically in the right panel.

3.4 Viscoelastic material properties

3.4.1 Rock

For rock, we used the viscosity and rigidity model used for exomoon interiors described in Dobos & Turner (2015) which is based on the works of Henning et al. (2009), Moore (2003), and Fischer & Spohn (1990). The behavior of the rock changes as a function of temperature, as the temperature crosses the solidus point (T_s = 1600 K), across a “breakdown” point (T_b = 1800 K) at which the volume fraction of solid rock crystals is equal to the volume fraction of melted rock. At this point, the behavior of the crystal and melt mixture begins to behave more like a liquid (Renner et al. 2000). When the rock is completely melted, (above the liquidus (T > T_l, where T_l = 2000 K) its viscosity and shear modulus are quite small. For temperatures below the solidus, μ_r = 50 GPa, and the viscosity is described by, $${\eta }_{\text{r}}(T)={\eta }_0\exp (\frac{Q_{\text{r}}^{*}}{R_{\text{G}}T}),$$(13)

where η₀ = 2.13 × 10⁶ Pa s (Henning et al. 2009), the activation energy $Q_{\text{r}}^{*}=333$ kJ mol⁻¹ (Fischer & Spohn 1990), and R_G = 8.314 J (mol K)⁻¹ is the gas constant.

For temperatures above the solidus, but below the breakdown temperature (T_s < T < T_b), the shear modulus depends on temperature (Fischer & Spohn 1990), $${\mu }_{\text{r}}(T)={10}^{\left(\frac{{\mu }_1}{T}+{\mu }_2\right)},$$(14)

where μ₁ = 8.2 × 10⁴ K and μ₂ = −40.6. In this regime, the viscosity of rock depends on the volume fraction of melt present (Moore 2003), f, which varies linearly with temperature between the liquidus and solidus, $${\eta }_{\text{r}}(T)={\eta }_0\exp (\frac{E}{R_{\text{G}}T})\exp (-Bf),$$(15)

where B is the melt fraction coefficient (10 ≤ B ≤ 40, Moore 2003). Above the breakdown temperature, but still below the liquidus, the shear modulus, now controlled by the behavior of the volumetrically dominant silicate melt, is a constant at μ_r = 10⁻⁷ Pa. The viscosity follows the Einstein-Roscoe relationship (Moore 2003), $${\eta }_{\text{r}}(T)={10}^{-7}\text{Pas}\exp (\frac{40000\text{K}}{T}){(1.35f-0.35)}^{-5/2}.$$(16)

Above the liquidus, μ_r = 10⁻⁷ Pa s, and the viscosity continues to depend on temperature, obeying Eq. (16) with f = 0.

3.4.2 Ice I and liquid water

For ice I, deformation was assumed to occur by volume diffusion (Goldsby & Kohlstedt 2001; Barr 2008), for which the viscosity depends on temperatureand grain size d, $${\eta }_{\text{I}}(T)=\frac{R_{\text{G}}T{d}^2}{14V_{\text{m}}{D}_{\text{ov}}}\exp (\frac{Q_{\text{ice}}^{*}}{R_{\text{G}}T}),$$(17)

where the molar volume V _m = 1.95 × 10⁵ m³ mol⁻¹, diffusion coefficient D_ov = 9.10 × 10⁻⁴, and activation energy $Q_{\text{ice}}^{*}=59.4$ kJ mol⁻¹. Grain sizes for ice in natural systems on Earth are limited to a few millimeters, due to the presence of silicate particles which inhibit grain growth (Durand et al. 2006). Temperature and strain rate conditions in terrestrial glaciers are somewhat similar to those expected in the outer planet satellites in our solar system, leading to the suggestion that d ≲ 1 mm in the ice I mantles of those bodies (Barr & McKinnon 2007). In the absence of other information, it is reasonable to assume that the process would work in a similar fashion on other planets: we adopted d = 1 mm as a nominal value. Similar to rock, the shear modulus for ice I is weakly temperature-dependent (Gammon et al. 1983), $${\mu }_{\text{I}}(T)={\mu }_{I,T_{\text{m}}}[\frac{1-aT}{1-a{T}_{\text{m}}}],$$(18)

where ${\mu }_{I,T_{\text{m}}}=3.39\times {10}^9$ Pa, a = 1.418 × 10⁻³ K⁻¹.

The melting temperature of ice I is T_m = 273 K. Above the melting point, μ_I = 0, and η_I = 10⁻³ Pa s which we used for the liquid water layers.

3.4.3 High-pressure ice polymorphs (HPPs)

Experimental data for the ductile behavior of the high-pressure ice polymorphs are more limited than experimental studies for ice I. It is known that each of the high-pressure polymorphs has a distinct flow law, and many of them are strongly non-Newtonian, so that the viscosity depends strongly on stress (Durham et al. 1997b). Similar behavior is observed in rock far from the melting point (e.g., Karato et al. 1995). We used a flow law that approximates the laboratory-determined flow behavior forices VI and VII from Durham et al. (1997a,b). Our tidal model assumed a Maxwell behavior, in which viscosity does not depend on stress, so we assumed a constant stress of 0.3 MPa. The approximate hpp flow law, $${\eta }_{\text{hpp}}=A\exp \left(\frac{Q_{\text{hpp}}^{*}}{R_{\text{G}}T}\right),$$(19)

where A = 6.894 × 10⁻⁷ and $Q_{\text{hpp}}^{*}=110$ kJ mol⁻¹, gives a viscosity at the approximate melting point of ices VI and VII, T_m = 285 K, of 10¹⁴ Pa s, somewhat similar to nominal near-melting-point viscosities for ice I. Data for the shear modulus of ices VI and VII are also limited, but at temperatures close to the melting point, μ_hpp ≈ 3.5 × 10¹⁰ Pa (Shimizu et al. 1996), where this value represents an average over the range for which experimental exist, which roughly corresponds to the pressure conditions expected within the TRAPPIST-1 planets’ interiors. We assumed μ_hpp was independent of temperature below the melting point. Above the melting point, high-pressure ices take on the properties of liquid water: μ_hpp = 0, and η_hpp = 10⁻³ Pa s.

3.4.4 Iron

The viscosity and rigidity of iron were assumed to be constant, and equal to the values deduced for the inner core of the Earth. The rigidity of iron, μ_Fe = 1.575 × 10¹¹ Pa, based on PREM (Dziewonski & Anderson 1981). The viscosity of iron in Earth’s inner core is thought to lie between 2–7 × 10¹⁴ Pa s (Koot & Dumberry 2011); we used η_Fe = 4.5 × 10¹⁴ Pa s.

3.5 Benchmark

To test the validity of our approach, we first used our methods to calculate the equilibrium interior temperature and tidal heat flux from Jupiter’s moon Io, for which the tidal heat flux has been measured by spacecraft to be ~ 1 W m⁻² (Spencer et al. 2000; Veeder et al. 2004). We used the following parameters for Io: mass M = 8.93 × 10²² kg, rock volume fraction ϕ_r = 0.95, iron fraction ϕ_fe = 0.05 (Schubert et al. 2004), orbital period P = 1.769 days, eccentricity e = 0.0041, and R = 1821 km. With these values, our model yielded a tidal heat flux F_tidal = 0.38 W m⁻², comparable in magnitude to that observed, and an equilibrium rock mantle temperature T_eq = 1671 K, above the rock solidus, and consistent with the inference of a partially molten mantle beneath Io’s surface (Khurana et al. 2011).

3.6 Uncertainties

Because the uncertainties in the measured planetary radius and mass are quite large (in some cases, they permit zero mass), it is important to provide limits to our calculations, as well. Here, we have reported a range of interior structures for the lowest and highest densities that can be estimated from the minimum and maximum radii and masses (see Table 1). In the cases of planets e, f, and h, we have modified some values in order to avoid the possibility of zero mass and the lack of a rocky core. For these planets, when determining the interior structure for the low-density case (and also for the other density values in the case of planet h), we assumed the rock to be Prinn–Fegley rock, as described in Sect. 2.2. These calculations with the lowest and highest densities effectively “bracket” the range of possible interior structures.

With these new structures, we recalculated the tidal fluxes, as well. Also, the uncertainties of the orbital eccentricities propagate into the tidal heating calculation. For this reason, we recalculated the tidal fluxes using the lowest and highest eccentricity values given in Table 1. Hence, we got nine different values for the tidal flux (one mean value for each density and the upper and lower limits provided by the minimum and maximum eccentricities).

When calculating the incident stellar flux, the largest uncertainty comes from the unknown albedo of the planets. For these calculations, we used the albedo value of the Earth, A_B = 0.3. Because the true albedos of the planets are not known, we repeated the calculations for 0.1 and 0.5 planetary albedos, too, and used the difference as error bars for the stellar (and consequently for the global) flux of the planets.

For the runaway greenhouse flux, the uncertainty came from the change of the surface gravity which depends on the mass and radius of the planet which are highly uncertain. We used the maximum and minimum values of stellar flux obtained from the minimum and maximum density calculations, and indicate them as uncertainties of the global flux of the planets. The uncertainties of the tidal flux were neglected here, because those are two orders of magnitude lower compared to the uncertainties in the stellar flux.

4 Results and discussion

Table 2 summarizes the volume fraction of each constituent material in each of the TRAPPIST-1 planets. Because of the large uncertainties in the measured mass and radius of the planets, we first calculated three end-member structures, corresponding to the highest-density planet, a structure for the mean-density planet, in which the planet has the best-fit mass and radius from Wang et al. (2017), and the low-density planet corresponding to the lowest permitted density. For planet h, the mean density value was very close to the minimum that could be obtained with assuming a low-density rock core, thus, we reported only the mean and maximum density cases.

Given the large error bars on mass and radius for each planet, it was not always possible to determine the bulk composition or the number of different compositional layers inside the planet – for planets b, c, d, and e, water-free structures were permitted if the planets had densities close to their upper limits. Likewise, each of these planets could be composed of rock and H₂ O, completely devoid of iron, if their true densities were close to their lower limits. For planets b, c, d, and e, we also computed the limit where the structure qualitatively changed, for example, one of the layers disappeared or a new one appeared, consistent with each planet’s possible values of mass and radius.

The values in Table 2 served as inputs in our calculations of heat generation by tides and loss by interior processes. Table 3 summarizes the values of equilibrium mantle temperature, tidal heat flux, and global energy flux at the top of the atmosphere from our geophysical models. In addition to the uncertainties provided by the densities (mass and radius values), the orbital eccentricity altered the tidal heating flux on the planet, too. Calculating with the three density and three eccentricity values for each planet, we obtained nine different results for the equilibrium temperatures and the tidal heating fluxes. (For planet h it resulted only in six different values.) For simplicity, we only presented the one case calculated with the mean density and mean eccentricity value, and the upper and lower uncertainties from among all other cases. The equilibrium values of mantle temperature obtained in our calculations, as well as the tidal heat flux, ${\overline{F}}_{\text{glob}}$, and threshold heat flux for a runaway greenhouse, F_RG, are summarized in Table 3.

4.1 Planet b

The top row in Fig. 3 illustrates our proposed interior structures for TRAPPIST-1b. The mass of planet b is comparable to the planet Venus, but with a radius of 1.086 ± 0.035 R_⊕, its mean density of 3405 kg m⁻³ is significantly lower than the mean densities of either Venus or Earth. Given our assumed densities for rock and iron, the highest permitted density for TRAPPIST-1b, ρ_max = 5041 kg m⁻³, permitted the planet to be composed of these two materials, without a significant reservoir of H₂ O (see left panel in Fig. 3). However, this represents an extreme end-member scenario. The absolute lower limit of mean density, ρ_min = 2038 kg m⁻³, permitted an interior structure with no iron, solely composed of rock and H₂ O (see right panel in Fig. 3). Given the current estimates of mass and density, the true structure of TRAPPIST-1b lies somewhere between these end-members (for example, Fig. 3, middle), with an iron core, rock mantle, mantle of high pressure polymorphs, and a surface layer of H₂ O. This is consistent with more detailed interior models that take into account compression of rock and metal at depth, which suggest that TRAPPIST-1b must contain more than 6% to 8% H₂ O (Unterborn et al. 2018), likely about 25% H₂O by mass (Wang et al. 2017). It seems likely that the planet does have some liquid water, even if it is just a relatively thin layer at the surface (Unterborn et al. 2018). Due to its proximity to the star, the effective surface temperature on TRAPPIST-1b, T_eff= 400 K. Thus, the planet must be in a runaway greenhouse state.

The top left panel in Fig. 4 illustrates how the interior structure and composition of TRAPPIST-1b varies as a functionof planetary mass and radius. The figure depicts the range of masses and radii for which TRAPPIST-1b is composed of ice, rock, and metal in a four-layered interior structure similar to the middle panel in the top row in Fig. 3. At smaller masses and larger radii, TRAPPIST-1b is composed of ice and rock; for larger masses and smaller radii, the planet is made of rock and iron. If the error bars on the mass for TRAPPIST-1b were reduced to ± 0.1 M_⊕, it would be possible to determine which of the three compositions best matched the planet’s true composition.

We found that the equilibrium mantle temperature for TRAPPIST-1b was above the rock solidus temperature, suggesting that the rock component of TRAPPIST-1b can be partially molten. If TRAPPIST-1b has a thick mantle of solid hpp ice, buoyant melt will rise to the top of the mantle, where it can erupt at the rock/ice boundary. This will drive melting at the base of the high-pressure ice polymorph mantle, driving ice melting and the rise of buoyant water through the hpp ice mantle (Barr et al. 2001).

The surface heat flux from tidal heat for TRAPPIST-1b, F_tidal = 2.68 W m⁻², is twice that of Jupiter’s moon Io (Spencer et al. 2000; Veeder et al. 2004), which is one of the most volcanically active planetary bodies in our solar system (Schubert et al. 2004). Thus, we expect TRAPPIST-1b to be volcanically active, but any activity would be concentrated at the ice/rock boundary, or at the bottom of a thick surface ocean. This raises the possibility of habitable regions inside TRAPPIST-1b in the vicinity of, for example, hydrothermal vents at the ice/rock boundary, as has been suggested for Jupiter’smoon Europa (Vance et al. 2007). However, the global energy flux on TRAPPIST-1b is about three times larger than the runaway greenhouse limit for the planet, indicating that the surface temperature on the planet may be too high to support life.

Fig. 3 Three end-member structures for each of the TRAPPIST-1 planets, assuming constant densities for ice I and liquid water (white) with ρI = 1000 kg m−3, high-pressure ice polymorphs (ρhpp = 1300 kg m−3, light gray), rock (ρr = 5000 kg m−3, dark gray) or low-density rock as indicated in Table 2 (ρr = 3300 kg m−3, dark gray, see Sect. 2.2), and iron (ρFe = 12 000 kg m−3, black), corresponding to the maximum density permitted by orbital solutions (left), the mean density (middle), and lowest density (right). The sizes of the planets and thicknesses of ice polymorph, rock, and iron layers are quantitatively accurate; but in some cases, the thickness of the outer water/ice I layer has been expanded slightly to show detail.

Fig. 4 Boundaries between ice/rock, ice/rock/iron, and rock/iron compositions as a function of planetary mass and radius for the four innermost TRAPPIST-1 planets. Star indicates the location of the best-fit mass and radius from Wang et al. (2017). For TRAPPIST-1b (top left), gray regions show pairs of mass and radius for which the planet has a four-layer structure, similar to the “mean density” case illustrated in the top row in Fig. 3. For TRAPPIST-1c (top right), the vast majority of pairs of mass and radius values indicate a rock/iron composition. For planets d and e (bottom left and right), an ice/rock composition is most likely; structures completely devoid of H2 O are indicated only if each planet has a mass close to the highest permitted value and a radius close to the minimum permitted value.

4.2 Planet c

TRAPPIST-1c has the highest mean density among the planets in the system, with an estimated $\overline{\rho }=7642$ kg m⁻³. The simplest permitted structure with such a high density is a large iron core, with a rock mantle, and no H₂ O layer at the surface (see middle panel in the second row in Fig. 3). This structure is consistent with other interior models, which show that the planet is roughly 50% iron by mass (Wang et al. 2017), and less than 6% to 8% H₂ O by mass (Unterborn et al. 2018). Figure 3 also illustrates the two extreme end-member interior structures that were possible given the large error bars on the mean density for TRAPPIST-1c. Taking the highest possible density for TRAPPIST-1c, the simplest interior structure was essentially a massive iron core with a thin layer of rock at the surface. The other extreme possibility was a rock-rich planet with no metal core and a mantle of H₂ O (see right panel in Fig. 3).

The top right panel in Fig. 4 illustrates the pairs of mass and radius for TRAPPIST-1c where a rock/iron composition is permitted, and the small region of parameter space where TRAPPIST-1c could be composed of rock and H₂ O. Further improvements to the mass measurement of TRAPPIST-1c could help distinguish between these models; if we knew the mass to within ± 0.5 M_⊕, we could determine whether TRAPPIST-1c might harbor a layer of H₂O near its surface. Consistent with others’ interior models, we found that H₂ O-bearing structures for TRAPPIST-1c are possible only if its density is close to the lower limit permitted by the error bars.

The equilibrium mantle temperature in the rock layer inside TRAPPIST-1c was T_eq = 1666 K, above the solidus temperature for rock, and a surface tidal heat flux of 1.3 W m⁻², comparable to Jupiter’s moon Io. Because it is unlikely that TRAPPIST-1c has a surface layer of H₂ O, it is possible that TRAPPIST-1c may have surface eruptions of silicate magma. The high heat flow may be large enough to be observed in infrared wavelengths. The heat flow on Io is known to be spatially and temporally variable as new volcanic vents form, or eruptions occur from previously dormant calderas (e.g., Veeder et al. 2012). The surface of Io also harbors numerous lava lakes whose surfaces periodically rupture and founder, resulting in momentary increases in the global heat flow (Veeder et al. 2012). Thus, as instrumentation improves, searching for excess thermal emission from the TRAPPIST-1 system may provide an avenue for the discovery of volcanism in the system.

In the mean density case, the iron core had a radius of ~ 4850 km which is larger than the iron core of the Earth (radius ~3500 km). This impliesthat planet c might have a strong magnetic field protecting the surface from stellar wind and flare erosions. However, due to the close-in orbit, the global energy flux is almost twice as large as the runaway greenhouse limit, which probably results in an uninhabitable surface.

4.3 Planet d

Planet d is much smaller than its inner two neighbors, with a mass M = 0.33 ± 0.15 M_⊕, and radius R = 0.772 ± 0.030 R_⊕. The third row in Fig. 3 illustrates the range of possible interior structures permitted by the mass and radius estimates. In the high-density limit, the planet could be composed solely of rock and iron, with an iron core of radius 2820 km. In the low-density limit, the planet could have no iron and be composed of rock and H₂ O (right panelin Fig. 3). If planet d has any significant water content, it would be concentrated near the surface of the planet and would form a liquid water ocean at the surface because T_eff > 273 K.

Figure 4 illustrates the range of possible compositions for TRAPPIST-1d. For the vast majority of masses and radii, a significant H₂ O component was required to match the observations. If the mass of the planet is close to its maximum value, and radius close to the minimum permitted value, structures with no H₂O can match the observations. Determining the mass of TRAPPIST-1d to ± 0.07 M_⊕ would shed light on whether TRAPPIST-1d has significant water content.

The equilibrium interior temperature in the rocky component of TRAPPIST-1d T_eq = 1618 K was extremely close to the solidus temperature for rock, indicating that the mantle could be partially molten, but certainly low-viscosity and convective. Because TRAPPIST-1d likely has a thick layer of H₂O on the surface, any melt extraction from the mantle would result in melting of the solid ice mantle. With an effective surface temperature T_eff = 288 K, TRAPPIST-1d is likely to be covered by a global water ocean. With a global heat flux from tidal heating of 160 mW m⁻², TRAPPIST-1d experiences a tidal heat flux twenty times higher than Earth’s mean heat flux, which provides an abundance of geothermal energy to drive chemistry in its surface ocean.

We only found an equilibrium temperature for the high-density case for high orbital eccentricities. In the cases of lower eccentricity, the tidal heating rate was not high enough to drive convection in the body, hence no equilibrium temperature can be found. In other words, this indicated that the planet may not experience significant tidal heating, and/or that a more sophisticated tidal model may be needed for these cases to determine the balance between tidal heating and convection. The error bars for the eccentricity of planet d permitted the possibility that the planet’s orbit is almost circular, implying insignificant tidal heating.

According to the calculations of Kopparapu et al. (2017), planet d is in the runaway greenhouse state. In contrast, we found that the global energy flux is lower than the runaway greenhouse limit, if the planetary albedo is ≳ 0.3. However, Kopparapu et al. (2017) uses a climate model for Earth-like bodies, but our model Pierrehumbert (2010) took into account the mass of the planet (through its surface gravity), which was significantly smaller than the Earth (0.33 ± 0.15 M_⊕). Because the planet was very close to the runaway greenhouse limit (if not exceeded it), liquid water on its surface was only possible if the planet was tidally locked to the star. Because of the synchronization, a high albedo cloud above the substellar point can reflect most incoming radiation (Yang et al. 2013; Kopparapu et al. 2016).

4.4 Planets e, f, and g

As a group, TRAPPIST-1e, f, and g have somewhat similar masses, and radii, yielding similar mean densities. Each of the three planets is less massive than Earth but with a radius comparable to that of the Earth, yielding mean densities in between those for ice and rock, similar to the mean densities of outer planet satellites and dwarf planets in our solar system, including Pluto (Stern et al. 2015), Enceladus (Porco et al. 2006), and Ganymede (Pappalardo et al. 2004).

The mass of planet e is quite uncertain, $M={0.24}_{-0.24}^{+0.56}{M}_{\oplus }$: the error bars encompass M = 0. As a result, the mean density of planet e is also not well-constrained: $\overline{\rho }={1713}_{-413}^{+4790}$ kg m⁻³ (Wang et al. 2017). In the context of our assumed densities for ice, rock, and iron, the minimum density permitted for planet e was ${\overline{\rho }}_{\text{min}}=1300$ kg m⁻³, the density of the high-pressure ice polymorphs. With this density, the minimum permitted mass for planet e is M = 0.206 M_⊕, with a radius of R = 0.957 R_⊕. Assuming alow-density core composed of Prinn–Fegley rock with ρ_r = 3300 kg m⁻³, TRAPPIST-1e could harbor a small rock core, even though its minimum density was comparable to ρ_hpp (see Table 1). The fourth row in Fig. 3 illustrates the end-member interior structures, ranging from a water-free planet to one composed of almost 100% H₂ O (right). With an effective surface temperature T_eff = 251 K (Gillon et al. 2017) below the melting point of water ice at low pressures (273 K), it is possible that TRAPPIST-1e has a solid ice surface. The relative sizes of the liquid water ocean and polymorph layers depended on the temperature in the interior of the planet, but the maximum thickness of the ice I layer can be constrained by the depth to the 209 MPa phase transition among ice I, liquid water, and ice III (Hobbs 1974).

The bottom right panel in Fig. 4 illustrates the pairs of mass and radius for which TRAPPIST-1e had an H₂ O-rich composition. If the mass of TRAPPIST-1e is close to the maximum permitted value, the planet could contain only rock and iron, but for the vast majority of mass/radius pairs, a thick ice mantle is predicted. To rule out a water-free structure, one would have to determine the mass of TRAPPIST-1e to within ~ 0.5 M_⊕.

Planets f and g have masses and radii similar to planet e, and thus, similar mean densities and interior structures. Both planets are less dense than TRAPPIST-1e, and so it seems likely that both of the planets have a significant amount of H₂ O. Similar to planet e, the mass of TRAPPIST-1f is not well-constrained, and the error bars permit the planet to have a low-density structure with a mean density equal to that of the high-pressure ice polymorphs. Because there is a thin layer of ice I or water at the surface, the low-density structure could have a very small Prinn–Fegley rock core.

For planets e and f, we found that their rock mantles would have equilibrium temperatures near the solidus of rock, similar to planet d. Also their global heat fluxes from tidal heating were ~ 160 to 180 mW m⁻² similar to planet d. Both planets f and g are far enough from their parent star to have T_eff < 273 K (Gillon et al. 2017), suggesting that their surfaces could be solid H₂O ice. Given that the tidal heat flow is about twenty times stronger than Earth, it seems likely that both planets could harbor liquid water oceans, perched between the surface layer of ice I and underlain by a mantle of hpp ices, similar to the jovian moons Ganymede and Callisto (Schubert et al. 2004).

We were not able to determine an equilibrium temperature for planet g, nor for the low-density case of planet f. This could be because tidal forces from the star are too small to significantly deform the rock cores of these planets, implying that tidal dissipation in the rock mantles of each planet is negligible. Tidal dissipation could generate heat in the ice mantles of the planet, but a more sophisticated model for k₂ would be needed to determine magnitude of the heat and its partitioning between the hpp ice mantle and ice I shell.

4.5 Planet h

The outermost planet, TRAPPIST-1h, has a mass of only M = 0.086 ± 0.084 M_⊕, and a radius of 0.715 R_⊕ (Wang et al. 2017). Given our assumed densities of rock and ice, TRAPPIST-1h, with a mean density $\overline{\rho }=1299$ kg m⁻³, should contain only H₂O, and no rock or metal.

Although it was permitted in the context of our assumptions, we thought it unlikely that TRAPPIST-1h would be composed of pure water ice. The vast majority of ice-rich bodies in the outer solar system contain some fraction of rock (one possible exception being Saturn’s small moon Tethys, which may not be primordial (Salmon & Canup 2017)). In our solar system, almost all of the outer planet satellites and dwarf planets contain some rock, and from a formation standpoint, it seems unlikely that TRAPPIST-1h is completely devoid of rock. For this reason we determined interior structures using a low-density mineral assemblage commonly used as a proxy for the rocky component of outer solar system satellites: Prinn–Fegley rock, for the planet, and a planetary mean density of 1300 kg m⁻³.

We found that a model TRAPPIST-1h with a density ${\overline{\rho }}_{\text{mean}}$ could harbor a core of rock 1030 km in radius, with the remainder of the planet composed of H₂ O in various phases (see right panel in the last row in Fig. 3). If TRAPPIST-1h had a density close to its maximum permitted value, ${\overline{\rho }}_{\max }=3149$ kg m⁻³, it could be composed of almost pure rock, with only ~100 km of H₂O near the surface (see left panel in the last row in Fig. 3). With its low effective surface temperature, T_eff= 169 K (Gillon et al. 2017), planet h undoubtedly has a solid ice surface. Similar to our null result for TRAPPIST-1g, we were notable to determine an equilibrium temperature dictated by the balance of tidal heat generation and transport inTRAPPIST-1h. Further knowledge of the composition and interior structure of TRAPPIST-1g and h, as well asa more sophisticated tidal model, will be required to assess the role of tidal heating in its interior evolution.

5 Implications for habitability

TRAPPIST-1b may contain a large amount of H₂O, and is likely volcanically active, but its activity is probably concentrated at the ice/rock boundary deep inside the planet. Due to its proximity to the star, the planet is very likely to be in a runaway greenhouse state (Kopparapu et al. 2017; Bolmont et al. 2017; Bourrier et al. 2017). TRAPPIST-1c is not likely to have much water or ice, and could have surface eruptions of silicate magma in the style of Io. Detection of excess thermal emission from the TRAPPIST-1 system could provide indirect evidence for surface volcanism on these bodies.

Using 1-dimensional and 3-dimensional climate models, Gillon et al. (2017) deduce that the surface temperatures of planets e, f, and g are suitable for harboring water oceans on their surface. However, using the runaway greenhouse limit of Pierrehumbert (2010), we found that planet d, too, may be habitable, if its albedo is ≳ 0.3. Planet d might be covered by a global water ocean that can provide a favorable environment for the appearance of life. This conclusion supports the finding of Vinson & Hansen (2017), who calculate the surface temperature of the TRAPPIST-1 planets due to stellar irradiation and tidal heating due to the circularization of the orbits, and found that planets d, e, and f might be habitable.

Planets e, f, g, and h are likely to have icy surfaces with liquid oceans underneath. According to Checlair et al. (2017), however, if the planets are tidally locked and have an active carbon cycle, then they should not have solid ice surfaces. The three-dimensional climate models of Turbet et al. (2018) show that TRAPPIST-1e is able to sustain its water content in liquid phase regardless of its atmospheric composition, provided that the planet is synchronously locked and abundant in water. In the atmospheres of planets f and g, a few bars of carbon dioxide is needed to facilitate the melting of ice on their surfaces. Most studies find that planet e is most likely to be habitable among the TRAPPIST-1 planets. For example, according to Wolf (2017), even a thin atmosphere could be enough to sustain temperatures of ~ 280 K locally around the substellar point, melting ~13% of the surface ice. Vinson & Hansen (2017) also find that including tidal heating to the energy budget, the surface temperature of planet e is around 280 K.

According to the tidal heating calculations of Papaloizou et al. (2017), planet e, and in most of the investigated cases, also planet f, could be both in the conservative habitable zone and in the tidal habitable zone. The tidal habitable zone indicates an intermediate level of tidal heating flux on the planetary surface (between 0.04 and 2 W m⁻²) which could drive plate tectonics required for the carbonate-silicate cycle that stabilizes the atmospheric temperature of the planet (Barnes et al. 2009). Papaloizou et al. (2017) also find that planet d orbits in the optimistic habitable zone, and its tidal heating is sufficient to be in the tidal habitable zone, too.

Bourrier et al. (2017) find that during the assumed 8 Gyr lifetime of the system, planets b through f may have lost, together, more than 20 Earth oceans worth of water due to hydrodynamic escape. However, the estimation is highly uncertain because of the wide possible ranges in planetary masses, and their values should be considered upper limits for hydrodynamical water loss. Magnetic interaction with the host star, however, can result in additional loss of water (Dong et al. 2018).

In any case, on planets b, e, f, g, and h, life might appear in the tidally heated (subsurface) ocean close to hydrothermal vents. Although our tidal heating model did not constrain equilibrium temperatures for planets g and h, they may still experience tidal heating. Our model might not be suitable for low heating rates: tidal heating can also be balanced by conduction alone if the planet’s mantle is too viscous to convect. We have also neglected planet–planet tidal interactions which may significantly contribute to the tidal heating in each planet, because they orbit in a close mean motion resonant chain. In addition, using tidal evolution models, Barnes (2017) find that planets in the habitable zone of late type M dwarfs are very likely to enter a 1:1 spin–orbit resonance in less than 1 Gyr. If the planet’s rotation is fast, then the orbital eccentricity is expected to grow slowly, which will lead to somewhat stronger tidal heating.

Lingam & Loeb (2018) investigate the probability of habitability of three TRAPPIST-1 planets, e, f, and g, orbiting inside the circumstellar habitable zone. By calculating the likelihood functions for habitability based on surface temperature and atmospheric loss of the planets (due to hydrodynamic escape and stellar winds), they find that for planet e, which gave the most promising results, the likelihood function is at least one to two orders of magnitude lower than in the case of Earth.

Any form of life is threatened however, by strong stellar activity. Analysis of K2 observations of TRAPPIST-1 showed frequent flares which present a serious threat for planetary atmospheres (Vida et al. 2017; Garraffo et al. 2017). Strong magnetic fields are required to protect the atmosphere from erosion, and life from high energy radiation.

6 Summary

The seven planets of the TRAPPIST-1 system have been the focus of intense observation and debate since the discovery of the system. Several of the planets are roughly Earth-sized, and their mean densities are not dissimilar from the terrestrial planets and large icy satellites in our solar system. Moreover, the effective surface temperatures of each of the planets are relatively modest, suggesting solid or water-covered surfaces.

Each of the seven TRAPPIST-1 planets are quite close to the parent star, as indicated by their orbital periods which range from 1.5 to about 20 days. The planets also occupy a mean-motion resonance (Gillon et al. 2017) and have non-zero orbital eccentricities. This implies that the planets experience tidal heating, which can be a significant energy source.

Here, we have constructed interior models for the planets using constant densities for ice, rock, and metal. The masses and radii of each of the planets are not well-constrained and have large uncertainties; for this reason, we applied our model using the mean density of each planet and also the extrema obtained from the error bars of observations. Within the range of estimated masses and radii for planets b, c, d, and e, each planet could be composed of solely iron and rock (our “high-density” cases). If the true densities of planets b, d, and e are close to the mean values implied by observations, each of these planets could have a layer of H₂O near the surface. The large density of planet c points toward a rock/iron composition, H₂ O is permitted only if the true density is closer to the lower limit of the estimated value.

Planets b and c likely harbor magma oceans in their interiors. The surface heat fluxes from tidal heat range from hundreds of milliwatts per meter squared (about twenty times higher than the Earth’s heat flow) to watts per meter squared, comparable to the heat flow estimated for Jupiter’s volcanically active moon, Io. Further work is required to reconcile the seemingly water-rich composition of planet b with its proximity to the TRAPPIST-1 star. The majority of permitted compositions for planet c indicate that the body is water-poor, raising the possibility that melt extracted from the magma ocean could erupt onto the surface. Further observations of the TRAPPIST-1 system may be able to detect the thermal or chemical signature of this activity, if present.

Of the seven planets in the TRAPPIST-1 system, we consider planets d and e to be the most likely candidates for habitable environments. Planet d is likely to be covered by a global water ocean, and can be habitable if its albedo is ≳0.3. Planetary albedos may be estimated from photometric measurements during occultations (see e.g., Rowe et al. 2008) which may give an additional constraint on the habitability of TRAPPIST-1d. As discussed in Sect. 5, planet e is likely to have liquid water on its surface, too, according to climate models (Turbet et al. 2018; Wolf 2017).

Planets f, g, and h are undoubtedly H₂O-rich, with thick mantles of ice, possible subsurface liquid water oceans, and solid surfaces composed of ice I, the low-density phase of ice present in terrestrial glaciers and ice sheets. It is worth noting that the present estimates of the mean density of planet h permit a 100% H₂O composition, but this seems unlikely given that even the most ice-rich bodies in our solar system (for example, the Kuiper belt objects) have a significant amount of rock in their interiors (Barr & Schwamb 2016). The geodynamical evolution of the outer planets of TRAPPIST-1 should be studied in more detail with a layered tidal and thermal model that can self-consistently determine the amount of tidal heating and heat transport by convection in both the rock and ice mantles.

Given the uncertainties in the compositions of each of the planets, we have employed a relatively simple geophysical model to determine the balance between tidal heat production and transport in each planet. Our approach provides more realistic values for Im (k₂) than a constant-Q approach because we take into account the feedback between the temperature of the planet. We used a tidal model that assumes each planet is a homogeneous sphere, but mimic the effect of the presence of multiple materials by determining an effective viscosity and rigidity for the planet based on theviscosities and rigidities of its component materials. We focused our efforts on determining the temperature of the rock mantle in each planet by balancing heat produced with the amount of heat transported by solid-state convection in the rock layer of each planet. When updated masses and radii are available, we will be pursuing more detailed calculations of tidal heating using multi-layered models and numerical techniques.

Further observations of the TRAPPIST-1 system are required to provide better constraints on the composition and interior structure of each planet. When improved compositional constraints are obtained, more sophisticated modeling techniques that calculate tidal dissipation and energy transport in a multi-layered planet can provide improved estimates of the tidal heat flux and mantle temperature. Our knowledge of the compositions of the four interior planets, chiefly, their water content, can significantly be improved if their masses can be determined to within ~ 0.1 to 0.5 M_⊕, which will require further observations of the TRAPPIST-1 system.

Acknowledgements

We thank the anonymous referee for a thorough review of the manuscript and helpful comments. We also thank Dr. László Molnár for initiating this collaboration, and Dr. Jessica Irving, Dr. Tess Caswell, and Dr. Ákos Kereszturi for helpful conversation. V.D. has been supported by the Hungarian National Research, Development and Innovation Office (NKFIH) grants K119993, K-115709 and GINOP-2.3.2-15-2016-00003.

References

All Tables

All Figures

Fig. 1 Schematic representation of our tidal and thermal model. Left: internal structure for the planet, composed of water ice I (white), high-pressure ice polymorphs (HPP, gray), rock (dark gray), and iron (black). The volume fractions of each constituent material were calculated based on the mean density of each planet (see Sect. 2). Middle: to calculate the response of each planet to tidal forcing, we constructed a volume-weighted average viscosity and rigidity for the planet based on the ϕ values and the rheology of each constituent material. A value of the imaginary part of the k2 Love number, and the tidal heat flux Ftidal were calculated. Right: the temperature in the planet’s interior (Teq) was governed by a balance between heat generation by tides and transport by solid-state convection; the efficiency of both processes depends sensitively on the viscosity and rigidity of the planet. The net global heat flux ${\overline{F}}_{\text{glob}}$ was calculated from the tidal heat flux and the incident sunlight from the star, F⋆.

In the text

Fig. 2 Left: schematic representation of the heat flux (in log units) from tidal heating (black line) and convective heat transport (gray) as a function of temperature. Two equilibria are possible: point (A) represents a balance between tidal heating and convection in a purely solid mantle (Teq is less than the solidus temperature Ts). This stateis depicted schematically in the middle panel. Point (B) represents an equilibrium for which Teq > Ts, implying a partially molten rock mantle exists within the planet. This state is shown schematically in the right panel.

In the text

Fig. 3 Three end-member structures for each of the TRAPPIST-1 planets, assuming constant densities for ice I and liquid water (white) with ρI = 1000 kg m−3, high-pressure ice polymorphs (ρhpp = 1300 kg m−3, light gray), rock (ρr = 5000 kg m−3, dark gray) or low-density rock as indicated in Table 2 (ρr = 3300 kg m−3, dark gray, see Sect. 2.2), and iron (ρFe = 12 000 kg m−3, black), corresponding to the maximum density permitted by orbital solutions (left), the mean density (middle), and lowest density (right). The sizes of the planets and thicknesses of ice polymorph, rock, and iron layers are quantitatively accurate; but in some cases, the thickness of the outer water/ice I layer has been expanded slightly to show detail.

In the text

Fig. 4 Boundaries between ice/rock, ice/rock/iron, and rock/iron compositions as a function of planetary mass and radius for the four innermost TRAPPIST-1 planets. Star indicates the location of the best-fit mass and radius from Wang et al. (2017). For TRAPPIST-1b (top left), gray regions show pairs of mass and radius for which the planet has a four-layer structure, similar to the “mean density” case illustrated in the top row in Fig. 3. For TRAPPIST-1c (top right), the vast majority of pairs of mass and radius values indicate a rock/iron composition. For planets d and e (bottom left and right), an ice/rock composition is most likely; structures completely devoid of H2 O are indicated only if each planet has a mass close to the highest permitted value and a radius close to the minimum permitted value.

In the text