© The Authors 2026

1 Introduction

Since the first detection of an exoplanet orbiting a solar-type star over two decades ago (Mayor & Queloz 1995), the number of detected exoplanets has increased rapidly, with several thousand now identified¹. While initial efforts focused on exoplanet detection, the field is now shifting towards their characterization. Although the detection of small temperate exoplanets has been challenging, about 40 have been identified to date with some characteristics similar to Earth (e.g., Turbet et al. 2023), including those in the TRAPPIST-1 system (Gillon et al. 2017). This system provides an excellent opportunity to study the atmospheres and surfaces of rocky planets outside our solar system (e.g., with JWST, Greene et al. 2023; Zieba et al. 2023), which could revolutionize our understanding of the evolution and habitability of terrestrial planets through comparative planetology (e.g.,Turbet et al. 2020).

The TRAPPIST-1 (also referred to as T1) system consists of seven Earth-sized planets orbiting a Jupiter-sized (very low mass) star forty light-years away (Gillon et al. 2017). At least three of these bodies orbit within the classical habitable zone of the star. These planets have semi-major axes ranging from 1 to 6% of an astronomical unit, forming a compact and stable planetary system due to their resonant architecture (Luger et al. 2017). An inherent particularity of the system due to its interacting planets is that, while transit photometry brought precise measurements of the planets’ sizes (Delrez et al. 2018; Ducrot et al. 2020), their strong mutual interactions made it possible to measure very precisely their masses, hence giving access to their densities. Their composition is consistent with rocky compositions, more volatile-rich or iron-poor than the Earth’s, at least for the outer planets (Grimm et al. 2018; Dorn et al. 2018; Agol et al. 2021; Park et al. 2025).

Internal heating in rocky bodies shapes their interior and surface characteristics as well as their evolution. Among internal heating sources, tidal dissipation can represent a large source of energy for planetary interiors. The most striking evidence in the Solar System is the case of Io, archetype of tidally-heated world hosting extreme volcanism (e.g., Spencer et al. 2000; Kervazo et al. 2021; Pettine et al. 2024; Mura et al. 2025). Interestingly, the orbital periods and eccentricities in the TRAPPIST-1 system are similar to those of the Galilean satellites, and both systems have mean-motion resonances. Tidal dissipation is therefore thought to contribute significantly to the total energy budget of the TRAPPIST-1 planets, with tidal heat fluxes being at least an order of magnitude larger than the Earth’s mean heat flux (Turbet et al. 2018; Bolmont et al. 2020a), and it is expected to lead to the persistence of magma oceans in the closest TRAPPIST-1 planets (Barr et al. 2018; Dobos et al. 2019). Given the recent claim of tidal heating occurring on the LP791-18d planet (Peterson et al. 2023), and the ongoing JWST observations of the TRAPPIST-1 system, it is now timely to reevaluate the estimate of the tidal heat flux for these planets.

This work thus generalizes to the work of Bolmont et al. (2020a), which was focused on the outer planets of TRAPPIST-1. Here, the use of the BurnMan code² to compute the internal structure of the TRAPPIST-1 planets allows us to investigate the tidal response of the inner planets and to investigate the impact of degeneracy on the internal structure of the planets. Furthermore, we use up-to-date estimations of the radius, mass, and eccentricities of the TRAPPIST-1 planets (Agol et al. 2021).

We present our models, both the internal structure model and the tidal model, in Sect. 2. We then give new estimates accounting for uncertainties in the internal structure and eccentricities of the planet in Sect. 3.

2 From internal structure to tidal heating

In this section we give the characteristics of both the internal structure model (Sect. 2.1) we used and the model that allows us to derive (from the internal structure model) the gravitational Love number of degree 2 and thus the dissipation in the planet (Sect. 2.2). Finally, Sect. 2.3 explains the method for computing the volumetric tidal heating inside the planets.

2.1 Internal structure model

For a solid planet of given mass and radius, interior models estimate the mass fractions of its building materials – typically refractory elements like iron and silicates, as well as volatiles such as water and organics – distributed among a metallic core, a rocky mantle, and an external envelope of varying mass and thickness. The composition of the planet has been proposed to be correlated with the metallicity of the host star (e.g., Unterborn et al. 2018), but it might not be a perfect correlation with the composition of rocky planets that span a wider distribution than that of stars (Plotnykov & Valencia 2020). In this work, we therefore adopt an agnostic strategy and consider a wide range of Fe/Si, resulting in a wide range of core sizes. As a result, in order to study a wide range of solid exoplanet interior compositions, several core sizes have been envisioned.

For each TRAPPIST-1 planet, we consider masses and radii from Agol et al. (2021). Despite the high precision of these estimations, internal structures and compositions are still degenerate and can be fit with various models. In this study, we assume a silicate mantle with an Earth-like composition and a liquid metallic core characterized by a Fe, Si, and S ratio. To account for the degeneracy, we change the relative size of the core of the planets compared to the mantle, and we calculate their composition accordingly in order to reproduce both their mass and radius. Our reference case is a planet that has the same core composition as Earth, characterized by a Fe/Si mass ratio of 1.2 (e.g., McDonough & Sun 1995). To keep parameter space manageable, we do not investigate volatile-rich interior compositions, although these could also reproduce the planets’ masses and radii (Agol et al. 2021).

We construct the internal structure models with the BurnMan code (Cottaar et al. 2014; Myhill et al. 2022), which allows us to integrate consistently physical properties such as density, temperature, gravity, and pressure from the surface to the planet’s center, assuming hydrostatic equilibrium. This tool lets us compute self-consistently depth-dependent density, temperature, shear modulus, and viscosity of each planet through a layered model consisting of, from the center to the surface, a liquid metallic core, a pyrolitic mantle, and a dunite rigid crust.

BurnMan computes the 1D profiles of thermoelastic and thermodynamic properties with the Birch-Murnaghan finite-strain equation of state (EoS). To do so, it uses lookup tables generated with Perple_X (Connolly 2005) to access information on properties such as heat capacity, conductivity, and thermal expansivity computed through thermodynamical equilibrium. BurnMan then uses this table to compute consistent profiles (density, thermal gradient, gravity gradient, P, and S wave, among others) for different pressure and temperature ranges. For a few cases we computed, the pressure at the bottom of the mantle was higher than the maximum pressure of the built-in lookup table included in BurnMan (135 GPa). In those cases, we expanded the lookup table using Perple_X. While BurnMan can technically extrapolate EoS values beyond experimental calibrations, we caution that uncertainties increase at high pressures, particularly for silicates, due to the lack of experimental constraints at high pressure. More details on the way internal temperature profiles are calculated can be found in Appendix A.

Viscosity is a key parameter in determining tidal heating efficiency in planetary interiors (e.g., Bolmont et al. 2020a). In our model, viscosity is computed from a classical Arrhenius law, temperature dependent: $$\eta =\frac{1}{2}{A}_0^{-1}{d}^{2.5}exp\left(\frac{E_a+P{V}_a}{RT}\right)$$(1) where V_a is the activation volume, E_a and A₀ are parameters depending on the material and d is the grain size. Equation (1) is valid for diffusion creep. In this study, the viscosity prefactor and grain size are kept fixed for simplicity, although they are poorly constrained and could vary by orders of magnitude. Here, we consider dry olivine: E_a = 300 kJ.mol⁻¹ and A₀= 6.08 × 10⁻¹⁹ Pa⁻¹.s⁻¹. We consider here a grain size d equal to 0.68 mm so that the factor $1/2A_0^{-1}{d}^{2.5}={10}^{10}\mathrm{Pa}.\mathrm{s}$. For the activation volume, we follow the best fit solution pressure dependent proposed by Stamenković et al. (2011). Using this prescription, typical values of V_a span a range of 1.8 cm³/mol (lower mantle) to 3.3 cm³/mol (upper mantle) for T-1b. As expected, this method leads to very high values of viscosity near the surface, so we set the viscosity of the rigid crust to 10²⁵ Pa.s. The viscosity of the inviscid liquid core is set to zero.

We first calculated a 1D “Earth-like” profile, which corresponds here to the internal structure the planets would have if they had the same composition as the Earth for the outer core (“Earth-like” column in Table 1), mantle, and lithosphere, from the center to the surface, respectively. The mantle was assumed to follow a perturbed adiabatic profile consistent with efficient internal cooling. This included a basal thermal boundary layer (TBL) with a contrast of 840 K. While plausible under Earth-like tectonics, such a TBL could be reduced in stagnant-lid planets or under strong tidal heating, which may suppress the gradient near the core-mantle boundary (CMB). The liquid iron core was assumed to follow an adiabatic profile and extend from the planet’s center to the core-mantle boundary. As we said previously, we prescribed the mantle as a pyrolitic mantle. Finally, the lithosphere was assumed to have a dunite composition with a thickness of 200 km.

In the reference model, we assumed a surface temperature of 300 K, a moho temperature (at the base of the crust) of 620 K, and a temperature at the base of the lithosphere (lab) of 1550 K. Then for each planet, we changed the amount of iron (Fe), silicon (Si), and sulfur (S) in the core, setting as a maximum amount of lighter elements, values constrained through chondritic models on Mercurian core (Vander Kaaden et al. 2020) and Mercurian models with eutectic composition (Harder & Schubert 2001). In this way, we considered the smallest core possible (which therefore had the highest iron content possible to be able to reproduce the mass of the planet), the biggest core possible (which therefore had the lowest iron content possible to be able to reproduce the mass of the planet), an intermediate case, and a Mercury-like case (where the core-size ratio was 85%). The extreme core case corresponds to a core size ratio even larger than the biggest core size measured in the solar system, which is why we also considered a Mercury-like case. Note that for the Earth’s core, Fe=90%, Si=2%, S=0% (e.g., Hirose et al. 2021). For the different compositions, we varied Fe, Si, and S while keeping the sum of Fe, Si, and S at 92%.

We also considered different surface temperatures for all planets: 300,600, and 800 K, and a temperature which should be more representative of the surface temperature for each planet (if they have an atmosphere made of CO₂ or O₂ for instance, following Lincowski et al. 2018). These reference temperatures are: 670 K for TRAPPIST-1b (or T-1b) and c, 650 K for T-1d, 300 K for T-1e, 250 K for T-1f, 210 K for T-1g, and 170 K for T-1h (M. Turbet, private communication, and compatible with Lincowski et al. 2018). The increment in surface temperature was then passed on to the moho and lab temperatures so that if the surface temperature was increased by 300 K, so were the moho and lab temperatures. For the different possible surface temperatures, we calculated the profiles following the procedure described above for the different cases (smallest core, intermediate core, Mercury-like case, biggest core), ensuring that we always reproduced the observed radius and mass. As the temperatures we considered were relatively similar, they had a small impact on the internal structure, so the core compositions we obtained were sometimes very similar for different surface temperatures. This can be seen in Table 1, where all the corresponding core compositions can be found for the different core size fractions and surface temperatures. The profiles of all planets assuming an Earth-like composition and their corresponding reference surface temperatures can be seen in Fig. B.1³.

Figure 1 shows the internal structure density profile for TRAPPIST-1b for the different temperatures and the core size assumptions. The major differences can be seen for the different core sizes. In particular, when the core is small, both shear modulus and viscosity have a much wider range in the mantle. The increase in shear modulus is due to the increase in pressure with increasing mantle thickness, while the slight increase in viscosity is related to the temperature and pressure profile computed with BurnMan. As a result, the viscosity is on average smaller if the planet has a big core than if it has a small core. Moreover, a bigger core and thinner mantle make the mantle more deformable.

While variations in core size significantly impact both shear modulus and viscosity, we find that changes in surface temperature have little impact on shear modulus but have a major impact on viscosity (top-right panel). Varying surface temperatures thus allow us to explore several viscosity distribution configurations, in accordance with Eq. (1). Due to this temperature dependency, there are almost three orders of magnitude between the viscosity we obtain at the base of the mantle between a surface temperature of 300 and 800 K. The different surface temperatures that we consider here, therefore, offer us a way to investigate the impact of different viscosity profiles on dissipation and the resulting tidal heating.

The rightmost panels of Fig. 1 display the silicate solidus temperature alongside the temperature profiles for TRAPPIST-1b (see Appendix B for details and applications to the other planets). We find that the temperature profiles for TRAPPIST-1b remain below the solidus throughout the mantle for the range of surface temperatures explored. However, profiles with a surface temperature of 800 K approach the solidus at depth, suggesting that even a modest increase in internal heating could trigger partial melting. While a full treatment of partial melting and its feedback on tidal dissipation is beyond the scope of this study, we note that approaching or crossing the solidus could substantially modify the planet’s tidal response. As shown by Kervazo et al. (2021), the presence of partial melt can enhance tidal dissipation, particularly for melt fractions near the percolation threshold. Our present estimates focus on solid interiors and should therefore be regarded as lower limits of the actual tidal heat production.

Fig. 1 Density, shear modulus, viscosity, and temperature profiles of TRAPPIST-1b computed with BurnMan code for different compositions (listed in Table 1). Top: influence of temperature on profile. Bottom: influence of composition on profile. The black curve shows the same composition in all panels (Earth-like composition and a surface temperature of 800 K). The gray curve in the top panels shows the curve corresponding to the reference temperature (here 670 K). The rightmost panels also show the temperature of the solidus (red dashed line).

2.2 Calculating the frequency dependence of the Love number

To calculate the distribution of tidal dissipation in a spherical multilayer body, we use the same method here as in Tobie et al. (2005), Bolmont et al. (2020a) and Kervazo et al. (2021). The method uses the elastic formulation of spheroidal oscillations developed by Takeushi & Saito (1972). Similarly to Dumoulin et al. (2017), we use the static formulation of Saito (1974) for the liquid core. This method was adapted to the viscoelastic case by Tobie et al. (2005), using the correspondence principle (Biot 1954). It was recently used to study multilayer solid planetary interiors (Tobie et al. 2019; Bolmont et al. 2020a; Kervazo et al. 2021). We refer the reader to these publications, in particular Appendix A in Kervazo et al. (2021), for more details about the method.

We make the assumption that the viscoelastic response of the planets follows an Andrade rheology (Andrade 1910; Castillo-Rogez et al. 2011)⁴. The Andrade rheological model is described by four parameters: the elastic shear modulus, μ, the shear viscosity, η, and two additional parameters, α and β, which describe the transient response between the purely elastic and viscous response. The α parameter determines the frequency dependence of the viscoelastic response. Following Tobie et al. (2019) and Bolmont et al. (2020a), we assume a fiducial value of α of 0.25, which is a typical value that reproduces the dissipation function of the Earth’s mantle over a wide range of frequencies. However, to account for the uncertainty on this parameter, we also consider the values of α=0.20 and α=0.30. For the β parameter, following (Castillo-Rogez et al. 2011), we assume that β ∼ μ^α−1 η^−α, which reproduces the existing data from mechanical tests on olivine minerals (Tan et al. 2001; Jackson et al. 2002). Note, however, that more recent studies point out that significant uncertainties regarding the appropriate β values remain, and β may possibly vary between 0.01 and 100 × μ^α−1 η^−α (Bierson 2024; Amorim & Gudkova 2024).

For each planet, and the associated tested configurations in terms of core size fraction and surface temperature (Table 1), we calculate the dependence of the Love number on the excitation frequency. The excitation frequency ω is a linear combination of the mean motion of the planet n and its spin Ω. Here, we consider the rotation of the TRAPPIST-1 planets to be synchronous, their eccentricities small, and their obliquity zero (e.g., Turbet et al. 2018). In that case, the main excitation frequency is ω=n. The main excitation frequency for each T1 planet is shown in Table 1.

Figure 2 shows the dependence of the imaginary part of the Love number Im(k₂) for T-1b. As seen in the top panel, increasing the surface temperature (and thus decreasing the viscosity) has a strong impact on the Love number. As shown in Bolmont et al. (2020a), decreasing viscosity shifts the maximum of Im(k₂) to higher frequencies, which leads here to an increase of the value of Im(k₂) at the excitation frequency of the planet (shown as a vertical dashed black line). In particular, between a surface temperature of 300 and 800 K, Im(k₂) increases by about a factor 4 (for α=0.25). The shaded areas in Fig. 2 show the dependence of the imaginary part of the Love number on α. This parameter does not change the values of Im(k₂) for the lower frequencies, but has an impact on the frequency range after the maximum of Im(k₂) (see also Fig. 2 from Bolmont et al. 2020a) where the slope is directly dependent on α (the slope is more pronounced for a high α). The differences grow the farther the frequency is from the frequency of the maximum. This means that for the highest surface temperature (orange curve), for which the maximum of Im(k₂) occurs closer to the excitation frequency, the impact of α on Im(k₂) is minimum. Consistently, for the lowest surface temperature (black curve), for which the maximum of Im(k₂) occurs farther to the excitation frequency, the impact of α on Im(k₂) is maximum.

The bottom panel of Fig. 2 shows the impact of increasing the size of the core is to shift the frequency of the maximum dissipation to higher frequencies. This is expected as increasing the size of the core decreases the average viscosity of the mantle (see Fig. 1) and decreasing the viscosity has this effect (see Fig. 2 from Bolmont et al. 2020a, and top panel of this figure). We can also see the impact of the presence and size of the core on the maximum dissipation (around a frequency of 10⁻¹²−10⁻¹¹ rad.s⁻¹). This maximum increases for smaller cores. This is also in agreement with the expected dependence of the imaginary part of the Love number on the shear modulus of the planet: the higher the shear modulus, the higher the maximum of the Love number (see Fig. 2 from Bolmont et al. 2020a). Indeed, Fig. 1 shows that the shear modulus is higher for smaller cores. As regards the top panel, due to the respective position of the peak compared to the excitation frequency, the impact of α on Im(k₂) is higher for the smallest core and lower for the largest core.

As discussed earlier, Fig. 2 also shows the excitation frequency of T-1b as vertical black dashed lines. This frequency is higher than the frequency of the maximum of the imaginary part of the Love number, which is also the case for all the other planets including planet h, which has the lowest excitation frequency. At that frequency, we can see that increasing the size of the core and the surface temperature increases the imaginary part of the Love number (bottom panel of Fig. 2). Consequently, we expect the highest tidal heating for the planets with the largest core and the highest surface temperatures. Note that for the two smaller core cases (green and black curves), the dissipation at the frequency of T-1b is very similar, so we expect similar tidal heating for these two internal structures.

Fig. 2 Frequency dependence of imaginary part of Love number for T-1b, for different compositions and temperatures (listed in Table 1). Top: influence of temperature. Bottom: influence of composition. The black curve shows the same composition in both panels (Earth-like composition and a surface temperature of 300 K). The gray curve in the top panel shows the curve corresponding to the reference temperature (here 670 K). The excitation frequency of T-1b is shown as the vertical black dashed line, which corresponds to its orbital frequency. The shaded region illustrates the dependence of the imaginary part on α (bracketed between 0.20 and 0.30). At the frequency of the planet, the dissipation for α=0.20 is higher than for α=0.30.

2.3 Calculating the volumetric tidal heating

We use the same method as in Bolmont et al. (2020a) to compute the volumetric tidal heating h_tide, which relies on Eq. (37) from Tobie et al. (2005) and is valid for synchronous, non-oblique planets on slightly eccentric orbits (e ≲ 0.05) $$h_{\text{tide}}\left(r\right)=-\frac{21}{10}{n}^5\frac{R_{\mathrm{p}}^4{e}^2}{r^2}{H}_{\mu }Im\tilde{\mu },$$(2) where R_p is the radius of the planet, r is the radius at which the volumetric tidal heating is estimated, n is the orbital frequency, and e is the eccentricity of the orbit; H_μ represents the radial sensitivity to the shear modulus μ. It depends on the radial structure of the planet and on y_i functions, which are associated with radial and tangential displacements, radial and tangential stresses and the gravitational potential (Takeushi & Saito 1972). We refer the reader to Tobie et al. (2005); Bolmont et al. (2020a) for a more in-depth explanation of the meaning of these y_i functions. The $Im\tilde{\mu }$ in Eq. (2) is the imaginary part of the complex shear modulus. As in Tobie et al. (2005), we assume that there is no bulk dissipation (even if it might be an important effect in the case of partially molten interiors, see Kervazo et al. 2021), and that all dissipation is associated with shear deformation. In this formalism, $\mathrm{ism},Im\tilde{\mu }$ contains all the information about the dissipation. An important remark we can make here is that the radial dependency of the tidal heating is not dependent on the orbit and rotation of the planet, but on the internal structure only.

Here we model only the dissipative power associated with the stellar eccentricity tides, and neglect other sources of dissipation due to obliquity tides, spin libration (Frouard & Efroimsky 2017) or planet-planet tides (Hay & Matsuyama 2019). Their contribution to the total dissipated power is expected to be small compared to that of the main eccentricity tides. The obliquity of the planets has been shown to be very small (<1 degree, see Turbet et al. 2018), and the amplitude of the spin librations was also estimated to be small (≲ 1 degree, see Revol et al. 2024). Previous studies suggest the possibility of large obliquity (Millholland et al. 2024) and large chaotic spin libration events (Vinson et al. 2019; Shakespeare & Steffen 2023; Chen et al. 2023). However, the high-obliquity stable state found by Millholland et al. (2024) requires unrealistically low dissipation for rocky material. In addition, the results of Vinson et al. (2019); Shakespeare & Steffen (2023); Chen et al. (2023) rely on the constant time lag (CTL) model, which is not well suited for studying the rotation of rocky planets (Makarov & Efroimsky 2013). A dedicated study is thus necessary to ensure that the heat from obliquity tides and spin librations is truly negligible, using an N-body code, for instance (Bolmont et al. 2020b; Revol et al. 2024). Finally, the heat generated by planet-planet tides should be lower than 2.5 × 10⁻² times the contribution of eccentricity tides (Hay & Matsuyama 2019). In any case, what we calculate here can be considered a lower estimate of the total heating, which should be slightly higher if we were to account for all these contributions.

For the eccentricities, we consider two extreme cases coming from Agol et al. (2021). First, we consider what could be the minimum eccentricities for the different planets, which are the forced eccentricities derived from the transit timing variations (TTVs) analysis performed in Agol et al. (2021). These forced eccentricities can be considered as the eccentricities the planets would have if they had had time to be tidally damped. Thus, these eccentricities arise solely from planet-planet interactions. Second, we consider what could be the maximum eccentricities for the different planets. To compute these maximum eccentricities, we compute the quadratic sum of the forced eccentricities and the value of the free eccentricities for the 95th percentile (which corresponds to a 2σ upper value). The values of these eccentricities are given in Table 2.

To compute the total tidal heating, we can either integrate h_tide from Eq. (2) over the planet as follows $$P_{\text{tide}}={\int \int \int }_V{h}_{\text{tide}}\mathrm{d}V=4\pi {\int }_{R_{\text{CMB}}}^{R_p}{h}_{\text{tide}}\left(r\right)r^2\mathrm{d}r,$$(3) where R_CMB is the radius of the core-mantle boundary, or use a global formula such as Eq. (2) of Tobie et al. (2005) $$P_{\text{tide}}=-\frac{21}{2}Im\left(k_2\right)\frac{{\left(n{R}_{\mathrm{p}}\right)}^5}{G}{e}^2.$$(4) The tidal heat flux is then ${\Phi }_{\text{tide}}=P_{\text{tide}}/\left(4\pi R_{\mathrm{p}}^2\right)$. We have systematically checked the agreement of the global heat power between the two formulations of Eqs. (3) and (4), and we get an error of less than 1.7% between them.

3 Tidal heating

We first discuss our results in the framework of the previous section (synchronous rotation, no obliquity, and a small eccentricity), but in the second subsection, we propose a convenient way to calculate the heating profile for a generic case.

3.1 For synchronous rotation, no obliquity, and a small eccentricity

Figure 3 illustrates the impact of (1) the four surface temperatures considered (T_surf=300, 600, 670, and 800 K, from panels a to d respectively, (2) the various sizes of core of Table 1 (colored lines), (3) the different eccentricities of Table 2 (full and dashed lines), and (4) the α parameter (extent of the colored areas) on the heating profile inside T1-b⁵.

Figure 3 shows that increasing the size of the core increases the tidal heating in the mantle, with the maximum dissipation occurring at the base of the mantle, as observed in Bolmont et al. (2020a). We can see a difference of tidal heating at the base of the mantle of a factor of 3−4 between the smallest core case (red) and the biggest core case (blue) for all surface temperatures. This is consistent with Fig. 2 in Sect. 2.2, where the difference in Love number at the excitation frequency was about this order of magnitude. Figure 3 also shows the influence of the assumption on the eccentricity of the planet (full or dashed lines). The maximum eccentricity of planet b is 32 times higher than the minimum eccentricity (see Table 2) and this leads to a difference in the tidal heating of three orders of magnitude (as shown by the e² factor in Eq. (2)). Finally, the extent of the colored areas between the full and dashed lines allows us to visualize the impact of the α parameter, which corresponds to an uncertainty of a factor of 2 in the tidal heating for a surface temperature of 300 K and the smallest core case. Consistent with what was discussed in Sect. 2.2, this uncertainty due to α decreases for bigger cores. As we observe in Fig. 2, the impact of α decreases with increasing surface temperature. Although not negligible, especially for the lowest temperatures (higher viscosities), the uncertainty on α is lower than the uncertainty on the surface temperature (viscosity) and on the internal structure and much lower than the uncertainty on the eccentricity. While this is true for planet b, this might be different for the other planets, for which the eccentricity is much better constrained than that of planet b.

A more convenient way to visualize the dependency of the tidal heating on the different parameters is to calculate the resulting total heat flux (calculated by Eq. (3) or Eq. (4)). Figure 4 shows the tidal heat flux for all planets of TRAPPIST-1, calculated with Eq. (3). Once again, we represent the influence of the eccentricity (blue for low eccentricities and green for high) and of the internal structure (colored areas delimited by triangles). Figure 4a shows the influence of the α parameter (different transparencies, with α=0.30 being the more transparent and α=0.20 being the less transparent). Figure 4b shows the influence of the surface temperature (different transparencies with the lowest temperatures, so highest viscosities, being the more transparent and the highest temperatures, so lowest viscosities, being the less transparent). The tidal heat flux is compared to the tidal heat flux of Io (red dashed line, e.g., Spencer et al. 2000; Lainey et al. 2009), the value of Earth’s geothermal heat flux (black dashed line Davies & Davies 2010), and the Earth’s tidal heat flux due to dissipation in the mantle (black dash-dotted line, computed from the profile in Fig. 3). Note that on Earth, most of the tidal dissipation occurs in the ocean (e.g., Egbert & Ray 2000; Egbert & Ray 2003). The values of the tidal heating of all planets can be found in Table 3. For this table, we give the values obtained for α=0.25 and the most likely surface temperatures (670 K for b and c, 650 K for d, 300 K for e, 250 K for f, 210 K for g, and 170 K for h), and the error bars encompass the rest of the uncertainties, which means that the highest value proposed is the one obtained that maximizes everything. So, the maximum values correspond to the largest core, α=0.10, and the highest temperatures, and the minimum values correspond to the smallest core, α=0.30, and the lowest temperatures. All combinations of parameters are given in the data accompanying this article⁶.

Figure 4 shows that for planet b (and planet c), the minimum and maximum eccentricities are quite different (more than one order of magnitude), which leads to a huge difference in the tidal heat flux. The minimum values of the tidal heat flux are between the value for the heating of the Earth and the tidal heating of Io, which is the most volcanic body of our Solar System, while the maximum values are well above the tidal heat flux of Io. Note that such high tidal fluxes are very likely to melt a large fraction of the rocky mantle, which would significantly impact the internal profile of the planet (both shear modulus and viscosity), and hence significantly modify the tidal response. In some circumstances, melt accumulation may result in the formation of mushy layers where tidal dissipation may be strongly enhanced due to localized reductions in viscosity and shear modulus (e.g., Kervazo et al. 2021) and possibly, in case of extreme melt production, to the formation of magma oceans where other modes of dissipation driven by gravito-inertial waves may develop, analogous to dissipation in the Earth’s ocean (e.g., Egbert & Ray 2003; Tyler et al. 2015; Farhat et al. 2022; Aygün & Čadek 2024; Farhat et al. 2025). While a full treatment of partial or total melting and its dynamical feedbacks is beyond the scope of this study, we acknowledge its potential impact on the tidal response. In our models, temperature profiles are computed independently of the predicted tidal heating rates. However, for the most irradiated planets such as TRAPPIST-1b and c, the high volumetric tidal heating could drive temperatures toward or beyond the silicate solidus, leading to partial melting in the mantle. As shown in Kervazo et al. (2021), partial melt can enhance tidal dissipation, especially when the melt fraction is close to the critical melt fraction. Above this limit, however, the material may behave more fluid-like, reducing the efficiency of viscoelastic dissipation.

Figure 4 also shows that uncertainty on the internal structure also leads here to significant uncertainty on the tidal heat flux for planets b and c, although to a lesser extent than uncertainty on their eccentricity. A factor of about 4−5 separates the flux of the smallest core case from the flux of the largest core case. For planets d, e, f, g, and h, the uncertainty on the eccentricity is much smaller. This means that for these planets, uncertainty on the tidal flux is dominated by uncertainty on the internal structure (i.e., the size of the core).

To focus on the impact of α (Fig. 4a) and viscosity via the surface temperature (Fig. 4b), Figure 4a shows that the uncertainty on α in the Andrade rheology leads to an uncertainty slightly lower than the uncertainty on the internal structure. The difference in global tidal heating between α=0.30 and α=0.20 is a factor of 3–4 while the difference between the smallest core case and the largest core case is a factor of 4–5. Figure 4b shows the impact of the chosen viscosity profile (via the surface temperature parameter). The darker shaded area corresponds to the highest temperature for all planets considered here (800 K). The light shaded area corresponds to a surface temperature of 300 K, and the hatched area corresponds to the most likely surface temperature. The highest temperatures (lowest viscosities) lead to the highest tidal heat fluxes, and the uncertainty this entails is a factor of 3−4 for the largest core case (most dissipative structures) and a factor of 4–6 for the smallest core case (least dissipative structure).

We can compare the tidal heat fluxes calculated with the top of atmosphere (TOA) fluxes for each planet from Ducrot et al. (2020). These TOA fluxes are shown as full black circles in Fig. 4. For most planets, the TOA fluxes are much higher than the tidal heat flux. However, assuming the highest eccentricity possible and the largest core configuration for planet b and the highest surface temperature (the configuration that maximizes the tidal heat flux for a given α, taken to be 0.25 for Fig. 4b), we obtain a tidal heat flux which is the same as the TOA flux. This means that a potential atmosphere of T1-b could be heated as much from the top as from the bottom, which should have repercussions on the atmosphere itself. In any case, the tidal heat flux of planet b, and to a lesser extent of planet c, could be quite high and potentially observable on JWST emission or phase curve data of the system.

Figure 4 shows that the tidal heat flux of T-1g and h should be lower than the heat flux of the Earth (black full triangle) and as low as the tidal heat flux of the Earth (black triangle; note that it corresponds only to the rocky part, so does not include the oceanic tide). Taking into account other types of heating, such as radioactive heating (e.g., Unterborn et al. 2022), induction heating (e.g., Kislyakova et al. 2017), or flare-induced heating (e.g., Grayver et al. 2022), these could potentially have a strong impact on the heat budget of the outer planets. However, the inner planets should have a flux that is dominated by tidal heat flux (planets b to d). Indeed, of the other sources of heating, it seems that flare-induced heating might be the most important in the context of TRAPPIST-1; Vissapragada et al. (2022) showed that it is comparable to the energy released by Earth’s radioactive elements today (so it should be comparable to the full black triangle value).

Assuming minimum eccentricities for the T1 planets, planets d and e have a similar or even higher flux than planet c. This is due to the fact that the forced eccentricity of planet c is about one order of magnitude lower than that of planet d and e (see Table 2). This difference in eccentricity, together with a higher value of the imaginary part of the Love number for their respective frequencies, compensates for the fact that planets d and e are farther from the star.

Interestingly, T1-e, which might be the most suitable to sustain an ocean of liquid water (e.g., Wolf 2017; Turbet et al. 2018, 2020), has a tidal heat flux that could be of the order of magnitude of the heat flux of the Earth and higher. This has strong implications for the possible habitability of the planet, as such heat fluxes could sustain volcanic activity and associated hydrothermal activities at the seafloor, similar to what has been shown on Jupiter’s moon Europa (Běhounková et al. 2021). It could also favor secondary outgassing, replenishing the atmosphere and surface in CO₂ and H₂O and thus favoring long-term habitability (e.g., Godolt et al. 2019).

Fig. 3 Volumetric tidal heating profile for T1-b for different surface temperatures (and thus viscosity profiles): (a) 300 K, (b) 600 K, (c) 670 K, (d) 800 K. The different colors represent the different structures listed in Table 1. The areas delimited by the full (dashed) lines correspond to the minimum (maximum) eccentricities given in Table 2. The extent of the areas represents the sensitivity of the profile to α, with the lower (left) limit corresponding to α=0.30 and the upper (right) limit corresponding to α=0.20. These profiles were obtained with Eq. (2). The tidal heating profile of the Earth is shown in a dashed black line as in Bolmont et al. (2020a). Additionally, we represent areas delimited by faint dotted lines. These profiles are compatible with JWST observational constraints on the nightside temperature of the T1-b (291 K at 2σ, 322 K at 3σ from Gillon et al. 2025), which are here hypothesized to be equal to a tidal temperature. The lower left limit thus corresponds to 291 K, and the upper right to 322 K.

Fig. 4 Total heat flux for all planets. Left: calculated for a surface temperature of 300 K (or a high viscosity). The transparency of the colored areas represents the dependency on α, with the more transparent (lower values of tidal heating) corresponding to α=0.30. Right: calculated for an α=0.25. The lighter shaded region (lower values of tidal heating) corresponds to a surface temperature of 300 K (high viscosities). The darker shaded region (higher values of tidal heating) corresponds to a surface temperature of 800 K (low viscosities). The hatched region corresponds to the reference temperatures (670 K for b and c, 650 K for d, 300 K for e, 250 K for f, 210 K for g, and 170 K for h). The colored area delimited by triangles represents the uncertainty we have on the internal structure, for a given assumption of the eccentricity (blue: minimum eccentricity, green: maximum eccentricity). These values are compared to the tidal heat flux of Io (full red triangle), Earth’s heat flux (full black triangle), and the Earth’s tidal heat flux (black triangle). The top of atmosphere (TOA) fluxes coming from Ducrot et al. (2020) are shown as full black circles. Finally, we show as red squares recent observational constraints from the JWST (Gillon et al. 2025) for the maximum nightside temperature of T1-b, which we assume is a tidal temperature (see Sect. 3.2).

3.2 For non-synchronous rotation, non-zero obliquity, and larger eccentricities

Equations (2) and (4) are valid for a specific set of orbital and rotational parameters. In particular, here we used expressions that are valid for synchronous rotation, zero obliquity, and small eccentricities (Tobie et al. 2005). No expressions exist in the literature for a generic formula for the tidal heating profiles in multi-layered planets, probably due to the fact that it is extremely challenging to derive them. We therefore propose here a workaround for this limitation (inspired by Tobie et al. 2005).

The solution stems from the fact that the shape of the tidal heating profile depends only on the structure, while its amplitude of course depends on the excitation frequency and orbital and rotational parameters of the system. This means that if we have a way of calculating the global tidal heat flux of a planet for a generic case, it is possible to shift a profile previously computed under more restricted conditions so that its integral is equal to the global tidal heat flux. This allows us to assess a generic heating profile.

This is what we have done in Fig. 3. We used observational constraints on the nightside temperature of TRAPPIST-1b obtained through the joint measurement of the phase curve of planets b and c (Gillon et al. 2025). Given the amplitude of the reconstructed phase curve of planet b, and previous studies (Greene et al. 2023), it seems very likely that it has no atmosphere. Whether the planet is tidally locked or not (Vinson et al. 2019), the nightside should be very cold. If it is tidally locked, it means the nightside is permanently in the dark. However, some studies have found that the T1 planets might not be perfectly tidally locked. For instance, Vinson et al. (2019) proposed that their rotation is chaotic, while Revol et al. (2024) proposed that there is a slow drift of the substellar point, which therefore leads to a day-night cycle. However, arguments based on the thermal inertia of rocks and the proposed value for the sidereal day (69 yr for T1-b, according to Revol et al. 2024), show that the day-night cycle is too long for the nightside to retain some heat from its passage on the dayside. We can therefore assume that a nonzero temperature from the nightside would come from internal heating. We consider here that the internal heating could be due to tides (therefore neglecting other sources of heating such as induction heating, Kislyakova et al. 2017).

The constraints in maximum nightside temperatures are T_2σ = 291 K at 2σ, T_3σ=322 K at 3σ from Gillon et al. (2025). From these temperatures, we can obtain the corresponding global fluxes with ${\Phi }_{2\sigma }={\sigma }_{\mathrm{SB}}{T}_{2\sigma }^4$ and ${\Phi }_{3\sigma }={\sigma }_{\mathrm{SB}}{T}_{3\sigma }^4$, where σ_SB is the Stefan-Boltzmann constant. Figure 4 shows the values of these fluxes (Φ_2σ=407 W/m², Φ_3σ=610 W/m²). They are relatively similar to the heat fluxes obtained previously for maximum eccentricity. However, we find some cases for which the tidal heat flux we obtain is higher than these observational constraints. This means that combinations of structure, viscosity profile (surface temperature), α, eccentricity leading to higher flux values than Φ_3σ can probably be rejected. If we consider that the maximum eccentricity allowed is the actual eccentricity of planet b (green color in Fig. 4), α=0.20, and a high viscosity (300 K of surface temperature, Fig. 4a), we can reject the hot (low viscosity) biggest core profile (upper green triangle facing downwards). If we consider that the maximum eccentricity allowed is the actual eccentricity of planet b, α=0.25 (Fig. 4b), we can reject (1) all structures for the highest surface temperature, meaning lower viscosities (darker shaded region), (2) all structures except for the smallest core case for a surface temperature of 650 K (hashed region). Only the high viscosity structures are compatible with the observations, whatever the core size (lighter shaded region).

From these global flux values, we can then calculate the corresponding volumetric heat profile. For each given internal structure, we start from a heat profile computed following the method described in the previous section, Sect. 3.1 (for instance, for the lowest eccentricity and α=0.25), and we compute the corresponding global heat flux Φ_ref. We then compute the ratio of Φ_2σ and Φ_3σ to Φ_ref and multiply the heat profile by this value.

Figure 3 shows the resulting profiles for both observational constraints as the shaded area delimited by a semi-transparent dotted line: the left dotted curve corresponding to 291 K and the right dotted curve corresponding to 322 K. Once again, the different colors are for the different internal structures, and by construction these profiles all amount to the same global flux.

Here, we considered observational constraints for the global flux, but the global flux can also be computed from the Love number and the orbital and rotational parameters of the planet. The total energy can be written as the sum of the rotational energy and orbital energy, such as $$E_{\mathrm{tot}}=E_{\mathrm{rot}}+E_{\mathrm{orb}}=\frac{1}{2}I{\Omega }_{\mathrm{p}}^2-\frac{\mathcal{G}{M}_{\mathrm{p}}{M}_{\star }}{2a}.$$(5)

The loss of total energy is then evaluated with the time derivative as $${\dot{E}}_{\mathrm{tot}}=I{\Omega }_{\mathrm{p}}\frac{\mathrm{d}{\Omega }_{\mathrm{p}}}{\mathrm{d}t}+\frac{G{M}_s{M}_p}{2a^2}\frac{\mathrm{d}a}{\mathrm{d}t}=-{\dot{E}}_{\mathrm{tide}}=-P_{\mathrm{tide}},$$(6) which corresponds to the amount of mechanical energy converted into tidal heat. To be consistent with the way the profiles are computed here, the derivatives $\frac{\mathrm{d}{\Omega }_{\mathrm{p}}}{\mathrm{d}t}$ and $\frac{\mathrm{d}a}{\mathrm{d}t}$ can be obtained following Boué & Efroimsky (2019). Using this method, one could compute the tidal heating profiles for any eccentricity, rotation, or obliquity, given an internal structure, a value of α, and a reference profile.

These steps are crucial for the next steps of this study, which would take into account the impact of tidal heating on the interior structure of the planet and on its tidal dynamical evolution. Indeed, the tidal heating profile could be computed at each timestep of the integration of the orbit and rotation of the planet (for instance, using ESPEM or SPIROID, e.g., Revol et al. 2023) and it could then be used to recompute a consistent interior structure. This interior structure could then be used to compute a new tidal Love number, which would allow us to compute the next timestep of the evolution of the semi-major axis, eccentricity, rotation, and obliquity of the planet. This, however, is out of the scope of the present study.

4 Conclusion

We here give new estimates of the tidal heat flux of the planets using the latest estimates of radii and masses (Delrez et al. 2018; Ducrot et al. 2020; Agol et al. 2021) and accounting for Andrade rheology (Bolmont et al. 2020a). We also propose a way to evaluate the tidal heating profile in planets that have parameters in rotation, eccentricity, and obliquity, which are not restricted to synchronous rotation, small eccentricities, and zero obliquity. The method we propose relies on the fact that the radial dependency of the profile only depends on the internal structure and the Andrade parameter α, and that it is possible to compute a global heat flux in a generic way. This is a first step to be able to one day study the retroaction of tidal heating on the internal structure of a multi-layered planet and on its corresponding tidal evolution. In particular, our models currently do not account for the thermomechanical feedback of tidal heating, such as temperature increases leading to partial melting. This omission may underestimate the actual heat flux in the most dissipative cases, especially for TRAPPIST-1b and c, where our predicted volumetric heating may locally approach or exceed the silicate solidus.

Concerning the tidal heating estimates we provide here, we find that it could be higher than that of Io for planets b and c, and that planets up to planet f could have tidal heat flux higher than, or of the same order of magnitude higher than, Earth’s heat flux. We also show that the uncertainty on the tidal heat flux is mainly due to the uncertainty on the eccentricity for planets b and c, and mainly due to the uncertainty on the internal structure (size of the core and viscosity profile) for the other planets. The uncertainties due to the internal structure are such that even with a very well observationally constrained system, the internal structure degeneracy is one major hurdle standing in the way of precise estimation of the tidal heating of planets. Here we allowed for a wide range of compositions (though restricting ourselves to rocky compositions and therefore neglecting the fact that a volatile or water rich interior composition has been suggested, Agol et al. 2021), but it might be possible to constrain this composition based on the composition of the star (e.g., Dorn et al. 2018; Unterborn et al. 2018). Agol et al. (2021) actually propose internal structures with core mass fractions (CMF) of about 21 p m 4 wt% for all planets, which is close to the CMF we obtain for the Earth case and the smallest core case. However, even assuming that the size of the core is known, the uncertainty on the viscosity profile (which we investigate via the surface temperature) is a factor of 4 (between the extremes we considered), and the uncertainty on the tidal parameter α is a factor of 2 (between the extremes we considered).

The existence of coreless terrestrial exoplanets has also been envisioned (e.g., Valencia et al. 2007; Fortney et al. 2007; Seager et al. 2007; Sotin et al. 2007; Elkins-Tanton & Seager 2008; Valencia et al. 2010). While no such configuration is observed among the planetary bodies of our Solar System, the discovery of exoplanets has broadened the range of possible interior compositions, motivating models of simple end-members. Our study, however, does not address this scenario. Indeed, in the case of the inner TRAPPIST planets, the hot interiors predicted by models would rule out this possibility, as these planets would have undergone melting early in their evolution, likely leading to the formation of a metallic core. Additionally, Huang & Dorn (2025) recently showed that oxygen partitioning rules out coreless TRAPPIST-1 planets.

JWST observations (secondary transit depth or phase curve) of planets b and c could allow constraints on the tidal heat flux of the planets, especially if those do not have an atmosphere or have very tenuous atmospheres, as hinted at by the first JWST observations of the inner planets (Greene et al. 2023; Zieba et al. 2023). Thanks to recent JWST observations (Gillon et al. 2025), we can reject low viscosity structures for T-1b if its eccentricity is high. Low viscosities lead to tidal heat fluxes higher than the constraints we have on the maximum nightside heat flux of the planet. Additionally, TRAPPIST-1e experiences a tidal heat flux that could be compatible with volcanism or plate tectonics, which makes it an even more interesting astrobiological target. Future work should thus focus on self-consistent models that integrate tidal heating, interior melting, and thermal evolution, in order to better assess the long-term geodynamic and observational signatures of intense dissipation in close-in rocky exoplanets.

Data availability

This work has made use of the BurnMan code, which is available on https://geodynamics.github.io/burnman/. All interior profiles used in this study can be found on https://zenodo.org/records/14884378.

Acknowledgements

The authors would like to thank the anonymous referee for their constructive comments. The authors would like to thank Lena Noack for the interesting discussions on viscosity, which led to some restructuring of the manuscript. E.B. and A.R. acknowledge the financial support of the SNSF (grant number: 200021_197176 and 200020_215760). This work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation under grants 51NF40_182901 and 51NF40_205606. M.K. and G.T. acknowledge the financial support of ERC PROMISES project (grant number: #101054470). The computations were performed at University of Geneva on the Baobab and Yggdrasil clusters. This research has made use of the Astrophysics Data System, funded by NASA under Cooperative Agreement 80NSSC25M7105.

References

Appendix A More details about BurnMan computations

The core is assumed to follow an adiabatic temperature gradient, consistent with a fully liquid metallic core. To compute this profile, we rely on BurnMan’s built-in thermodynamic treatment of iron alloys. Specifically, we use the Fe-S EoS from Saxena & Eriksson (2015), which is appropriate for planetary interior conditions up to several hundred GPa. To approximate the effect of light elements (e.g., S, Si), BurnMan applies a correction to the molar volume, which proportionally reduces the density and indirectly modifies thermodynamic properties such as the thermal expansivity and heat capacity. These quantities are then derived self-consistently from the EoS.

For the mantle, the temperature profile is constructed by combining a perturbed adiabatic gradient in the convecting region with fixed thermal boundary layers at the surface and at the core-mantle boundary (CMB). The mantle itself is modeled using an Earth-like peridotitic composition, and its thermodynamic properties are obtained from BurnMan’s implementation of the third-order Birch-Murnaghan EoS, based on Perple_X-generated mineralogy and phase equilibria. The thermal profile assumes a total contrast of 900 K, with 60 K across the upper thermal boundary layer and 840 K across the basal boundary layer. This approach follows previous studies modeling Earth-like convecting planets and is consistent with a regime of vigorous mantle convection characterized by a Rayleigh number of 10⁷. Therefore, the CMB temperature and boundary layer thickness is not directly specified but results from the imposed thermal structure.

Appendix B Profiles for all planets assuming an Earth-like composition

Figure B.1 shows the different profiles of density, temperature and pressure of all the T1 planets assuming an Earth-like composition for the core. We compare the temperature profiles computed with BurnMan to the silicate solidus, following Hirschmann (2000) for pressures between 0 and 10 GPa, Herzberg et al. (2000) between 10 and 22.5 GPa, and Monteux et al. (2016) from 22.5 to 136 GPa. For the surface temperatures considered here, all planets except TRAPPIST-1d exhibit mantle temperatures below the solidus, indicating no partial melting under the assumed conditions. For TRAPPIST-1d, however, the temperature profile approaches or slightly exceeds the solidus near the base of the lithosphere, suggesting that limited partial melting could occur in this region. Further work would be required to assess how such localized melting could affect the planet’s tidal response and internal heat transport (e.g., Kervazo et al. 2021; Běhounková et al. 2021).

Fig. B.1 Profiles of density, temperature and pressure of the TRAPPIST-1 planets computed with the BurnMan code, assuming an Earth-like composition for the core and mantle. The solidus is shown in dashed lines.

Appendix C Tidal temperature

Figure C.1 shows the tidal temperature for each planet, calculated as (Φ_tide/σ_SB)^1/4, where Φ_tide is the tidal heat flux of Fig 4. The tidal temperature is compared with the equilibrium temperature as estimated in Ducrot et al. (2020) and the brightness temperature for planets b and c (Ducrot et al. 2020). We also compare for planet b the tidal temperatures we compute to the observational constraints on the maximum nightside temperature obtained by Gillon et al. (2025). The details can be found in Sect. 3.2.

Tides contribute to the temperature of the planets from a few Kelvin (for the outer planets) to more than a hundred Kelvin for the inner planets. Especially for T-1b, the maximum tidal temperatures are close to the maximum nightside temperature obtained from JWST constraints, and this can be translated into constraints on the degree of synchronization of the planets as well as on their obliquity (see Gillon et al. 2025).

Fig. C.1 Same as Fig. 4 but for tidal temperature. These values are compared with the equilibrium temperature of all planets (full black circles) and the brightness temperature of T-1b and c (full red circles) estimated in Ducrot et al. (2020). We also give the recently acquired brightness temperatures measured by the JWST: the temperature of planet b comes from Greene et al. (2023), and that of planet c from Zieba et al. (2023).

All Tables

All Figures

Fig. 1 Density, shear modulus, viscosity, and temperature profiles of TRAPPIST-1b computed with BurnMan code for different compositions (listed in Table 1). Top: influence of temperature on profile. Bottom: influence of composition on profile. The black curve shows the same composition in all panels (Earth-like composition and a surface temperature of 800 K). The gray curve in the top panels shows the curve corresponding to the reference temperature (here 670 K). The rightmost panels also show the temperature of the solidus (red dashed line).

In the text

Fig. 2 Frequency dependence of imaginary part of Love number for T-1b, for different compositions and temperatures (listed in Table 1). Top: influence of temperature. Bottom: influence of composition. The black curve shows the same composition in both panels (Earth-like composition and a surface temperature of 300 K). The gray curve in the top panel shows the curve corresponding to the reference temperature (here 670 K). The excitation frequency of T-1b is shown as the vertical black dashed line, which corresponds to its orbital frequency. The shaded region illustrates the dependence of the imaginary part on α (bracketed between 0.20 and 0.30). At the frequency of the planet, the dissipation for α=0.20 is higher than for α=0.30.

In the text

Fig. 3 Volumetric tidal heating profile for T1-b for different surface temperatures (and thus viscosity profiles): (a) 300 K, (b) 600 K, (c) 670 K, (d) 800 K. The different colors represent the different structures listed in Table 1. The areas delimited by the full (dashed) lines correspond to the minimum (maximum) eccentricities given in Table 2. The extent of the areas represents the sensitivity of the profile to α, with the lower (left) limit corresponding to α=0.30 and the upper (right) limit corresponding to α=0.20. These profiles were obtained with Eq. (2). The tidal heating profile of the Earth is shown in a dashed black line as in Bolmont et al. (2020a). Additionally, we represent areas delimited by faint dotted lines. These profiles are compatible with JWST observational constraints on the nightside temperature of the T1-b (291 K at 2σ, 322 K at 3σ from Gillon et al. 2025), which are here hypothesized to be equal to a tidal temperature. The lower left limit thus corresponds to 291 K, and the upper right to 322 K.

In the text

Fig. 4 Total heat flux for all planets. Left: calculated for a surface temperature of 300 K (or a high viscosity). The transparency of the colored areas represents the dependency on α, with the more transparent (lower values of tidal heating) corresponding to α=0.30. Right: calculated for an α=0.25. The lighter shaded region (lower values of tidal heating) corresponds to a surface temperature of 300 K (high viscosities). The darker shaded region (higher values of tidal heating) corresponds to a surface temperature of 800 K (low viscosities). The hatched region corresponds to the reference temperatures (670 K for b and c, 650 K for d, 300 K for e, 250 K for f, 210 K for g, and 170 K for h). The colored area delimited by triangles represents the uncertainty we have on the internal structure, for a given assumption of the eccentricity (blue: minimum eccentricity, green: maximum eccentricity). These values are compared to the tidal heat flux of Io (full red triangle), Earth’s heat flux (full black triangle), and the Earth’s tidal heat flux (black triangle). The top of atmosphere (TOA) fluxes coming from Ducrot et al. (2020) are shown as full black circles. Finally, we show as red squares recent observational constraints from the JWST (Gillon et al. 2025) for the maximum nightside temperature of T1-b, which we assume is a tidal temperature (see Sect. 3.2).

In the text

	Fig. B.1 Profiles of density, temperature and pressure of the TRAPPIST-1 planets computed with the BurnMan code, assuming an Earth-like composition for the core and mantle. The solidus is shown in dashed lines.
In the text

Fig. C.1 Same as Fig. 4 but for tidal temperature. These values are compared with the equilibrium temperature of all planets (full black circles) and the brightness temperature of T-1b and c (full red circles) estimated in Ducrot et al. (2020). We also give the recently acquired brightness temperatures measured by the JWST: the temperature of planet b comes from Greene et al. (2023), and that of planet c from Zieba et al. (2023).

In the text