© G. Tobie et al. 2019

1 Introduction

The number of detected exoplanets with radius and mass comparable to that of the Earth is increasing drastically (e.g. Bonfils et al. 2013; Batalha et al. 2013; Marcy et al. 2014; Dressing & Charbonneau 2015; Coughlin et al. 2016; Dittmann et al. 2017; Gillon et al. 2017). Most of the low-mass planet candidates orbit at relatively close distances from their central star and are therefore likely subjected to large tidal forces due to the gravitational tide raised by the star. Planets with masses lower than 10 M_⊕ and with short orbital periods (<10–20 days) seem especially abundant around M-dwarf stars (e.g. Bonfils et al. 2013; Dressing & Charbonneau 2015). For example, the TRAPPIST-1 system recently discovered by Gillon et al. (2016, 2017) exhibits several small planets orbiting a low-mass star at relatively close distance (<0.1 AU), corresponding to orbital periods of a few Earth days. For such planets, tidal friction which controls both the rotational and thermal evolution, is expected to have a large impact on climate stability and the planetary heat budget, and hence on planet habitability (e.g. Lammer et al. 2009; Barnes 2017; Turbet et al. 2018). Therefore, properly quantifying the way in which the planets respond to tidal forces is crucial to assess their evolution and address their potential habitability.

The impact of tidal interactions on planetary evolution depends on the orbital configuration of the system, but also on the way planetary bodies deform under the action of tides raised by their central star, and in some cases by a potential planetary companion, like in the Earth-moon system. In most studies, the tidal response is described by using two simple parameters, the potential Love number, k₂ and the dissipation function, Q⁻¹. The potential Love number corresponds to the ratio between the variation of potential induced by tidal deformation of the interior and the external potential raised by the central star, while Q⁻¹ represents the fraction of energy which is dissipated anelastically. In many models studying tidal evolution of planetary rotations and orbits (e.g., Mardling 2007; Bolmont et al. 2014; Barnes 2017), global dissipation is classically parametrized using prescribed constant values of k₂ and Q⁻¹ or the time lag Δt between the onset of tidal forcing and the response of the body (with Δt = arcsin(Q⁻¹)∕χ, χ being the tidal angular frequency). However, the tidal response depends on the planet composition, structure and thermal state as well as on the frequency of the tidal forcing. This requires taking into account the specificity of interior models and of orbital configurations when predicting the appropriate tidal parameters and their time evolution.

Many models have been developed to compute the tidal response of a variety of planetary objects of the solar system in the past using multi-layer methods (e.g. Alterman et al. 1959; Kaula 1964; Sabadini et al. 1982; Segatz et al. 1988; Tobie et al. 2005; Wahr et al. 2009; Rivoldini et al. 2011; Beuthe 2013; Dumoulin et al. 2017). Most studies dedicated to exoplanets have used simplified approaches to predict the tidal response assuming, for instance, the formula derived for homogenous viscoelastic interiors (Henning et al. 2009; Efroimsky 2012; Makarov & Efroimsky 2014; Barr et al. 2018; Makarov et al. 2018), and thus neglecting the rheological effects of radial heterogeneity induced by the increase of temperatureand pressure with depth and the coupling between the different internal layers. Taking the latter parameters into account is particularly important at the extreme pressures and temperatures relevant to exoplanets, and when very different materials (metals, silicates, water) in both solid and liquid states are considered (e.g. Takeushi & Saito 1972; Saito 1974; Henning & Hurford 2014; Dumoulin et al. 2017). To our knowledge, only the study of Henning & Hurford (2014) used a multi-layer method to predict the tidal response of a variety of exoplanets. However, these latter authors considered interior models with constant values of density, rigidity and viscosity for each internal layer. Here we propose to compute the tidal response of exoplanets for a variety of compositions taking into account realistic interior models, including variations of mechanical properties with depth.

With the growing number of detected exoplanets, a variety of structural models have been developed this last decade to predict the internal structure of massive rocky planets and large icy worlds (e.g., Valencia et al. 2006; Sotin et al. 2007; Seager et al. 2007; Grasset et al. 2009; Swift et al. 2012; Weiss & Marcy 2014; Dorn et al. 2015). These studies showed thatthe mass–radius relationship is sensitive to planet composition, meaning that accurate determination of planet mass and radius combined with stellar elemental abundances may potentially provide constraints on the core size, mantle composition, and water fraction of the planet. Even if the uncertainties are still relatively high on the mass determination of small-size planets (e.g., Weiss & Marcy 2014), we can already observe some compositional tendencies. Rogers (2015) noticed, for instance, that most planets with a radius larger than 1.6 R_⊕ have density lower than terrestrial planets suggesting a significant fraction of volatile material, in the form of H/He atmosphere and/or water layers. Based on 2025 Kepler planets, Fulton et al. (2017) noticed that planets tend to prefer radii of either ~ 1.3 R_⊕ or ~ 2.4 R_⊕, with relatively few planets having radii of 1.5–2.0 R_⊕, suggesting two distinct planet classes: super-Earths below this radius transition and mini-Neptunes above it. More recently, Wu (2019) showed that this radius transition shifts systematically to a larger radius for planets around more massive stars. Most of the detected planets have so far been found around low-mass stars. While taking into account the stellar mass dependence, Wu (2019) showed that the masses of super-Earth planets would peak around 8 M_⊕ for planets around Sun-like stars, suggesting that super-Earth planets may extend up to 10 M_⊕. All the calculations and derived scaling laws that we provide in the present study apply for solid (rocky or icy) planets with no extended external envelop. Therefore, for planet radii between 1.5 and 2 R_⊕ and masses ranging between ~6 and 8 M_⊕, our proposed scaling laws should be used with care, after verifying that the derived density and estimated surface temperature of the planets are consistent with solid rocky or ice-rich planets.

In the present study, we perform a systematic computation of the tidal response for planet masses ranging between 0.1 and 10 M_⊕, with variable iron and water ice content. For that purpose, we use a numerical code initially developed for computing the tidal response of icy moons (Tobie et al. 2005), recently applied to Venus and Earth (Dumoulin et al. 2017). In Sect. 2, we determine the internal structure of a planet as a function of its mass and composition following the method of Sotin et al. (2007). The rheological properties are adapted from the Earth case assuming an Andrade viscoelastic rheology and are extrapolated to high-pressure ranges (see Sect. 3). The method for computing tidal deformation is presented in Sect. 4. Validation of the structural and tidal computation is then presented for the Earth case in Sect. 5. Section 6 displays the results for rocky planets, with various Fe/Si ratios, while Sect. 7 is devoted to ice-rich planets. For both types of planets, we show that the tidal Love number can be described as a function of planet mass using a relatively simple relationship. We also present the sensitivity of the dissipation function to planet mass, tidal frequency and assumed viscoelastic parameters. Some applications of our derived scaling laws to a selection of exoplanets are finally provided in Sect. 8.

Fig. 1 Computed interior structure for a given mass and bulk composition. The bulk composition is defined relative to a reference composition corresponding to the solar value for the Fe/Si and Mg/Si ratio and to the terrestrial value for the Mg/(Mg+Fe) ratio in the silicate mantle. For rocky planets, the water ice mass ratio is 0%, while it can reach values up to 50% for ice-rich planets. For the examples shown here, an ice mass ratio of 20% is considered and planet masses of 1 and 10 M⊕ are displayed. The interior structure is divided into three main layers: ice mantle, rock mantle and iron core. The dashed lines indicate the separation between the low-pressure and high-pressure phases for both ice and silicate mantle.

2 Computation of the interior structure

We consider two classes of planets: rocky Earth-like planets, and planets enriched in water ice. For the two planet classes, we consider planet masses ranging between 0.1 and 10 M_⊕. For a given mass and assumed composition, the surface radius as well as the radius of each internal interface is computed using the method described below. The interior is assumed to be differentiated into a metallic core (assumed fully liquid), a silicate mantle, divided in an upper and lower mantle, and, for ice-rich planets, a thick ice layer (see Fig. 1). We do not consider the possible presence of external fluid envelops (dense atmosphere, magma ocean or liquid water ocean) and leave these particular cases to be addressed in future studies.

The interior structure is modeled using the approach of Sotin et al. (2007). The same equations of state (EoS) and parameters are used for all layers, except for the surface liquid water layer which is replaced by low-pressure ices. The density profile in the upper part of the silicate mantle (P < 25 GPa) as well as in the low-pressureice layers (P < 2.2 GPa) is computed with a Birch-Murnaghan EoS, up to third order in finite strain. For simplicity, no phase transition is considered in the upper rock mantle (P < 25 GPa) and the low-pressure ice layers. A Mie-Grunëisen-Debye EoS is employed for the liquid iron core, the lower part of the silicate mantle, and the high-pressure ice layer (Ice VII). The iron-rich core is assumed to be entirely liquid, and therefore no inner solid core is considered here. In all layers, the temperature profile is assumed to follow an isentropic temperature profile.

The planetary composition is described in terms of bulk ratios of Fe/Si and Mg/Si, the Mg content of the silicate mantle, and the mass fraction of H₂O. For simplicity, we consider a Mg content in the silicate mantle equal to that of the Earth (defined as the mole fraction Mg/(Mg+Fe) = 0.9) for the two classes of interior models. The reference values for Fe/Si and Mg/Si are fixed to the solar values,0.977 and 1.072, respectively. For the rocky planets, the iron fraction is determined by the quantity: ${\delta }_{\text{Fe}}=\frac{(\text{Fe/Si})\text{-}{(\text{Fe/Si})}_{\oplus }}{{(\text{Fe/Si})}_{\oplus }}$, which varies between −50 and +50%. For ice-rich planets, the iron content ratio is fixed to the reference value (δ_Fe = 0%) and only the ice fraction (x_ice = M_ice∕M_tot) is varied between 0 and 50%.

Figure 2 displays the mass–radius relationship we obtained by computing a range of interior models with various iron and ice contents and the comparison with scaling laws using the polynomial coefficients provided in Table 1.

Fig. 2 Mass–radius relationship, normalized to the Earth, obtained for rocky and ice-rich planets with various iron contents (−50% ≤ δFe ≤ + 50%) and ice fractions (0 ≤ xice ≤ 50%), relative to the composition of Earth. Each point corresponds to the results obtained using the modeling approach described above, based on Sotin et al. (2007). The lines correspond to the scaling laws using the polynomial coefficients provided in Table 1.

3 Rheological parameters and assumptions

The viscoelastic response of the interior to tidal forces is controlled by the density and the mechanical properties of each layer composing the interior, namely the elastic bulk, shear modulii, and viscosity. To model the viscoelastic response, we use an Andrade rheology, which is consistent with existing experimental data on silicate minerals and ices (Castillo-Rogez et al. 2011; Efroimsky 2012). Compared to the standard Maxwell rheology, the Andrade rheology provides an accurate description of the frequency dependency of both elastic and dissipative responses on a wide range of forcing frequencies and periods. The Maxwell rheology provides a reasonable description of viscoelastic behavior for forcing periods close to the Maxwell time (defined as the ratio between viscosity and the shear modulus), but strongly underestimates viscous dissipation (or overestimates Q factor) for forcing periods much smaller than the Maxwell time (e.g. Sotin et al. 2009). The mantle of the Earth, for instance, has a Maxwell time of the order of 10 000 yr, while the tidal forcing is of the order of days. Applying a Maxwell rheology for the computation of viscoelastic tides on Earth results in an overestimation of the Q factor by several orders of magnitude (see, e.g. Fig. 2 of Sotin et al. 2009), while, as we show in Sect. 5, an Andrade rheology with appropriate parameters is able to reproduce the observed Q value for the interior of the Earth, at diurnal tidal frequencies as well as at lower forcing frequencies.

Complex quantities taking into account this frequency dependency are thus defined from the elastic “unrelaxed” shear modulus μ, the elastic “unrelaxed” bulk modulus K, and the Newtonian viscosity. The complex compliance, which corresponds to the inverse of the complex shear modulus (μ^c (χ) = 1∕J(χ)), is given for the Andrade model by: $$J(\chi )=\frac{1}{\mu }-\frac{i}{\eta \chi }+\beta {\left(i\chi \right)}^{-\alpha }\Gamma \left(1+\alpha \right)\text{},$$(1)

where χ is the tidal frequency, α and β are parameters describing the frequency dependence and the amplitude of the transient response, and Γ is the gamma function. Comparison with available experimental data for rock and applications to the Earth (see Sect. 5) indicate that the Andrade model is a good approximation to describe the anelastic attenuation at tidal frequencies (Castillo-Rogez et al. 2011). For the α parameter, we explore a range of values varying between 0.2 and 0.3, which frames the typical value required to explain the Q factor of the mantle of Earth (see Sect. 5). For the β parameter, following the approximation of Castillo-Rogez et al. (2011), we assume that β ≃ μ^α−1η^−α, which is justified for olivine minerals (Tan et al. 2001; Jackson et al. 2002).

The elastic bulk isentropic modulus, K, is derived from the density and pressure profiles: $K=\rho \frac{\text{d}P}{\text{d}\rho }$. The shear modulus, μ, in each solid layer is then estimated from the bulk modulus and the pressure. For the silicate part, the following relationship is used: ${\left(\frac{\mu }{K}\right)}_{\text{sil}}=0.52-0.5\frac{P}{K}$ for P < 25 GPa, ${\left(\frac{\mu }{K}\right)}_{\text{sil}}=0.63-0.885\frac{P}{K}$ for P > 25 GPa, which accurately reproduces the shear modulus profile of the mantle of Earth (Dziewonski & Anderson 1981; Stacey & Davis 2008). For pressures above 130 GPa, we assume a similar relationship. For the ice layers, a similar relationship constrained from existing experimental data at high pressures (Polian & Grimsditch 1983) is used: ${\left(\frac{\mu }{K}\right)}_{\text{ice}}=0.6-0.9\frac{P}{K}$. We note that $\left(\frac{\mu }{K}\right)=0.6$ corresponds to a Poisson body.

The viscosity is assumed to increase due to phase transition (ice I > high-pressure ice; low-pressure silicate mantle > high-pressure silicate mantle > very-high-pressure silicate mantle) or due to composition change (ice > rock). For sake of simplicity, the viscosity is assumed constant in each sublayer. The viscosity in the ice layers typically varies within the range 10¹⁴ – 10¹⁶ Pa s for low-pressure ice (ice I to Ice VI) and 10¹⁶–10¹⁸ Pa s, and possibly even larger values (e.g. Durham et al. 1997), for high-pressure ices (ice VII). Based on estimates for Earth, the viscosity in the silicate mantle typically varies within 10¹⁹ – 10²¹ Pa s in the upper mantle and 10²²–10²³ Pas in the lower mantle (e.g. Čížková et al. 2012, and references therein). For larger planets, the silicate viscosity is expected to strongly increase above ~120 GPa, possibly up to 10²⁵−10²⁶ Pa s and more (e.g. Karato 2011; Tackley et al. 2013).

The behaviour of the viscosity at the very high pressure and temperature conditions of extra-solar planets is a current matter of debate (e.g. Mocquet et al. 2014). Stamenković et al. (2012) used a semi-empirical homologous temperature approach to scale enthalpy variations with melting temperature and to assess the very high pressure and temperature dependence of the viscosity. They obtained a factor of two reduction for the activation volume of Mg-perovskite between 25 GPa and 1.1 TPa, resulting in a viscosity increase by approximately 15 orders of magnitude through the adiabatic mantle of a 10-M_⊕ planet. Their results contrast with the conclusions of Karato (2011) that viscosity should decrease with pressure when the latter exceeds approximately 0.1 TPa. Tackley et al. (2013) argued that thermal convection through a self-regulation effect controls the viscosity values reached at high pressure and should be of the order of 10²⁵ – 10²⁶ Pa s and remain relatively constant with depth. Moreover, as shown in the case of Venus, the tidal response for interiors with viscosity that varies with depth is comparable to that computed for interiors with constant viscosity (equal to the average of the variable profile; Dumoulin et al. 2017), meaning that depth-variable viscosity can be reasonably neglected at first order.

For each layer, we define a reference viscosity which is equal to 5 × 10²⁰ Pa s in the low-pressure silicate mantle layer, 5.10²² Pa s in the high-pressure silicate mantle layer (corresponding to the upper mantle and lower mantle of Earth, respectively), and 5 × 10²⁵ Pa s in the very-high-pressure silicate mantle, existing for M ≥ 1 M_⊕ (see Table 2).

4 Computation of the viscoelastic tidal response

The viscoelastic deformation of the planet interior under the action of periodic tidal forces is computed following the method of Tobie et al. (2005), adapted to terrestrial planets in Dumoulin et al. (2017). The iron core is supposed to be fully liquid, inviscid, and incompressible, following the static formulation of Saito (1974), and all the other layers consist of viscoelastic compressible solids. Using the density profile and rheological properties determined using the methodology described in Sect. 3, the Poisson equation and the equation of motions are solved for small perturbations in the frequency domain for compressible media (assuming finite values for the bulk modulus, K) and assuming an Andrade viscoelastic rheology (e.g. Castillo-Rogez et al. 2011).

The potential Love number, k₂, characterizing the potential perturbation, and the dissipation function, Q⁻¹, are computed by integrating the six complex radial functions associated with the radial and tangential displacements, the radial and tangential stresses, the gravitational potential and the potential continuity, as defined by Takeushi & Saito (1972). The formulation of the spheroidal deformation developed by Takeushi & Saito (1972) was initially derived for the elastic case. Nevertheless, as explained in Tobie et al. (2005), the same formulation can be used for the viscoelastic case by solving it in the frequency domain and by defining complex shear and bulk modulii, equivalent to the elastic modulii used in the elastic problem. For more details, please see Tobie et al. (2005) where applications to icy moons using a Maxwell rheology are presented. Here we use the Andrade rheological model, which is more appropriate and accurately reproduces the dissipation function of the interior of the Earth as well as its frequency dependence over a wide range of frequencies (see Sect. 5).

For the fluid core, the simplified formulation of Saito (1974) relying on two radial functions, is employed assuming a quasi-static and nondissipative fluid material. The solution in the solid part of the interior is expressed as the linear combination of three independent solutions, which reduces to two independent solutions in the fluid part. The system of six differential equations is solved by integrating the three independent solutions using a fifth-order Runge-Kutta method with adaptive step-size control from the center (r = 0 km) to the planet surface (r = R_P) and by applying the appropriate boundary conditions (see for more details Tobie et al. 2005). For a given frequency, the complex Love numbers, $k_2^{*}$, is determined from the radial function of the potential at the planet surface (r = R_P), and the global dissipation function, Q⁻¹, is equal to the ratio between the imaginary part and the modulus of the complex Love number $k_2^{*}:Q^{-1}=\Im (k_2^{*})/\Vert k_2^{*}\Vert .$

5 Model validation on the Earth case

Figure 3 compares our synthetic profiles of density, ρ, bulk modulus, K, and shear modulus, μ, derived following the method described in Sect. 2 and using the reference parameters (Fe/Si = 0.977, Mg/Si = 1.072, Mg/(Mg+Fe) = 0.9), to the Preliminary Reference Earth Model (PREM; Dziewonski & Anderson 1981). The density and bulk modulus is accurately reproduced throughout the entire silicate mantle, with only some slight departure in the upper mantle as we neglect any phase transition occurring in this layer. Obviously, as we neglect the presence of a solid inner core, the density is underestimated in the inner part of the core, which has only minor effects on the tidal response.

For viscosities in the ranges of 10²⁰–10²¹ Pa s and 10²²– 10²³ Pa s in the upper and lower mantles, respectively, α parameters between 0.2 and 0.3, and using the PREM profiles for ρ, μ and K, we obtain a potential Love number for a semi-diurnal tidal period between 0.302 and 0.306, which is consistent with the observed value estimated between 0.304 and 0.312 (Table 3, Ray et al. 2001). At the period of 18.6 yr, we also obtained a k₂ value consistent with the value inferred from observations. Using the synthetic profiles, we obtain a slightly smaller value (0.299–0.302 at 12.42 h and 0.312–0.336 at 18.6 yr) due to the modest difference in the density profile related to the absence of a solid innercore in our synthetic model. This remains, however, very close to the value for Earth, the discrepancy being less than 1%.

As illustrated in Fig. 4, the modeled Q value is very sensitive to the assumed values of mantle viscosity and α parameter. At the period of semi-diurnaltide (M₂), for a lower mantle viscosity in the range 10²²–10²³ Pa s, consistent with existing geophysical constraints (Čížková et al. 2012), the Q value of Earth of about 230–360 inferred from satellite tracking and altimetry (Ray et al. 2001) can be reproduced for α values between 0.22 and 0.29. For this range of α, we verified also that we reasonably reproduce the frequency dependence of the Earth Q (see Table 3), inferred from the fortnightly M_f tide (13.66 days; Ray & Egbert 2012), the Chandler Wooble (433 days; Furuya & Chao 1996; Benjamin et al. 2006) and the 18.6-yr tide (Benjamin et al. 2006).

Fig. 3 Comparison between our synthetic profiles (solid line) of density, ρ, bulk modulus, K, and shear modulus, μ, derived following the method described in Sect. 2 with those of the Preliminary Reference Earth Model (PREM; Dziewonski & Anderson 1981) (dotted line), assuming an interior composition corresponding to the reference case (see Table 2).

Fig. 4 Computed Q factor of the solid mantle of Earth for the 12.42-h tidal period using the synthetic profile displayed in Fig. 3 for various values of viscosity for the lower high-pressure (HP) mantle, ηHP (with ηLP∕ηHP = 100) and of α parameter for the Andrade rheological model. The bold black line indicates the estimated value of Q for the interior of Earth and the gray area shows the uncertainty range estimated by Ray et al. (2001).

6 Results for rocky planets

Figure 5 displays the computed k₂ as a function of planet mass for three iron contents ranging between −50 and +50%, relative to the reference value. Planets with higher iron content can be seen to be characterized by higher k₂, which is explained by the larger size of the liquid iron core, resulting in a larger deformation of the silicate mantle. For all tested values of iron content, we find that the dependency of k₂ with mass can be scaled to the surface gravity relative to the Earth using the following relationship: $$k_2({\delta }_{\text{Fe}})=k_2^{1M_{\oplus }}({\delta }_{\text{Fe}})-\left[\frac{1-{\left(\frac{g}{g_{\oplus }}\right)}^2}{1+{\left(\frac{g}{g_{\oplus }}\right)}^2}\right]\times \Delta k_2({\delta }_{\text{Fe}})$$(2)

with g being gravity. The factor $k_2^{1M_{\oplus }}$ corresponds to the k₂ value for a mass equal to that of the Earth. Δk₂ corresponds to the correction that should be considered relative to the mass of Earth to reproduce the k₂ value for mass ranges between 0.1 and 10 M_⊕. For masses below 0.1 M_⊕ and above 10 M_⊕, our interior model is no longer valid as the relationship between μ, K, and ρ likely diverge from that based on the Earth. The scaling law we derived here for k₂ is therefore only valid in the range of 0.1–10 M_⊕.

We also derived the variations of the coefficients $k_2^{1M_{\oplus }}({\delta }_{\text{Fe}})$ and Δk₂(δ_Fe) as a function of iron content using a polynomial fit to degree 2 with x = δ_Fe. The coefficients derived from this polynomial fit are provided in Table 4. As shown in Fig. 5, this simple relationship reproduces k₂ relatively well over the whole range of masses tested here (0.1–10 M_⊕). The tidal frequency and assumed viscosity have only a very minor effect on k₂ (less than 1%) as the tidal response is dominated by the elastic response for the range of tidal period tested here (0.5–100 days).

For the dissipation factor, Q, we also observe a systematic increase as a function of planet mass as well as a dependency with iron content (which controls the ratio between mantle thickness and core radius). However, Q is much more sensitive to tidal frequency and mantle viscosity than the mass and internal structure (Fig. 6). As expected (e.g. Efroimsky 2012), Q decreases with increasing forcing periods and decreasing viscosity. For high viscosity values in the very-high-pressure rock layer, Q increases with planet mass as the proportion of the rock mantle with P > 130 GPa and hence the average viscosity of the mantle increases with mass. For planet masses below 1 M_⊕, we note also a strong decrease of Q, which is explained by a decrease of the silicate mantle viscosity as the low-pressure phase becomes more predominant with decreasing mass.

We did not find any straightforward scaling for the Q variations as a function of planet mass. Nevertheless, for a given planet mass, we found that the variation of Q as a function of tidal frequency, viscosity, iron content and α parameter can be reproduced with the following relationship: $$Q={10}^{(A+B\alpha )}\times {\left(\frac{\chi \eta }{{\chi }_{\text{ref}}{\eta }_{\text{ref}}}\right)}^{\alpha },$$(3)

with η_ref being the reference value for the viscosity profile (see Sect. 2) and χ_ref = 2π∕1day. As illustrated in Fig. 7, this scaling reproduces the Q factor very well for a wide range of viscosity, frequency and iron content. This scaling remains valid as long as the tidal period is much smaller than the Maxwell time, τ_M, defined as the ratio between the viscosity, η, and the shear modulus, μ. For the range of viscosity tested here, the Maxwell time is always larger than 100 yr (up to 1 million years for massive highly viscous exoplanets), which is much longer than the tidal periods we considered (<100 days).

Fig. 5 Computed k2 for rocky planets with mass ranging between 0.1 and 10 Earth’s mass and iron content varying between −50 and +50% relative to the reference value, for a tidal period of 1 day. The lines represent the fit based on Eq. (2) and the polynomial coefficient provided in Table 4.

Fig. 6 Computed Q factor for rocky planets with mass ranging between 0.1 and 10 Earth masses (a) for a tidal period of 1 day and three iron content values varying between −50 and +50% relative to the reference value, and (b) for the iron content of the Earth and three tidal periods of 1 (square), 10 (triangle) and 100 (diamond) days. Two extreme values are considered for the very-high-pressure rock viscosity (ηVHP): 5 × 1022 Pa s (gray) assuming no viscosity increase relative to the reference Earth lower mantle and 5.1025 Pa.s (black) assuming a factor 1000 increase relative to the reference, i.e. the lower mantle of the Earth. Curve symbols are the same as in Fig. 2.

7 Results for ice-rich planets

For ice-rich planets, we observed similar tendencies as a function of mass to those for rocky planets (Fig. 8). The dependencies of the k₂ coefficients (Eq. (2)) as a function of ice content can be well reproduced using a polynomial fit up to degree 3 with x = x_ice, corresponding to the mass fraction of ice relative to the total planet mass. The corresponding coefficients are provided in Table 4.

As shown on Fig. 9, Q varies only moderately with planet mass. It is mostly sensitive to the ice fraction, viscosity and forcing frequency. Like for the rocky planets, we noticed a systematic dependency with frequency and viscosity. However, for ice-rich planets, this dependency is more complex as the Maxwell time of the ice layer, τ_M = η∕μ, is much shorter and hence comparable to the tidal period. For ηχ values ranging between 10⁹ and 10¹⁰ Pa, Q reaches a minimal value, which varies between 2 and 7 depending on the ice content and planet mass (Figs. 10, 11). For ηχ ≫ 10¹⁰ Pa, Q follows a behavior comparable to the rocky planets with $Q~{(\eta \chi )}^{\alpha }$. When approaching values of typically 10¹⁰–10¹¹ Pa, the frequency–viscosity dependency changes to Q ~ ηχ. Furthermore, after the minimal Q value is reached, ηχ ≪ 10¹⁰–10¹¹ Pa, the ηχ-dependency becomes negative, $Q~{(\eta \chi )}^{-1}$, consistent with the prediction for a Maxwell model (e.g. Efroimsky 2012).

This indicates a progressive change from an Andrade-like behavior to a Maxwell behavior with decreasing ηχ. For forcing periods between a few days and a few tens of days and viscosity between 10¹⁶ and 10¹⁸ Pa s, the values of ηχ typically ranges between 10¹⁰ and 10¹² Pa s, corresponding to an intermediate regime between Andrade and Maxwell behaviors. A purely Maxwell regime can be reached only for very long forcing periods (>100 days) or low viscosity values (≤10¹⁴–10¹⁵ Pa s).

We find that the Q factor variations as a function of ηχ can be qualitatively reproduced by the following relationship: $$Q=\frac{Q_{\text{min}}/2}{X}+{\left(\frac{2}{Q_{\text{min}}X}+\frac{1}{Q_{\alpha }}\right)}^{-1},$$(4)

where Q_min is the minimal value of Q, $X=(\eta \chi )/{(\eta \chi )}_{Q=Q_{\text{min}}}$ and Q_α is a function depending on the α-parameters similar to Eq. (3). However, as the function Q_α for ice-rich planets varies as a function of planet mass and composition in a more complex manner than for rocky planets, we failed to derive a good fit as a function of ice content comparable to what was obtained for the rocky planet (see Fig. 7). As displayed in Fig. 11, for rocky planets, the ηχ values are in a range where the frequency–viscosity dependence is controlled by the α parameter of the Andrade model (Eq. (1)), while ice-rich planets are in a regime in transition between a Maxwell and Andrade behavior, making the frequency–viscosity dependence more complex. For ice-rich planets, Q is lower than 100 with minimal values as low as 2, while for rocky planets Q typically varies between 100 and 1000, depending mostly on the planet mass and internal viscosity (Fig. 11).

Fig. 7 Scaling of the Q factor as a function of mantle viscosity η, tidal frequency χ, Andrade α parameter, and iron content, δFe, relative to reference values (ηref = 5 × 1022 Pa s, χref = 2π∕1day), for four planet masses (1, 2, 5 and 10 M⊕).

Fig. 8 Computed k2 for ice-rich planets with masses ranging between 0.1 and 10 Earth masses and H2 0 mass fraction up to 50% relative to total mass. The lines represent the fit based on Eq. (2) and the polynomial coefficient provided in Table 4.

Fig. 9 Computed Q factor for ice-rich planets with masses ranging between 0.1 and 10 Earth masses for two tidal periods (1 day in black and 100 days in gray) and three values of ice fraction (10% – cross, 20% – triangle, and 50% – diamond) for a HP ice viscosity of 1018 Pa s and thereference viscosity profile for the rocky part (see Sect. 2). Curve symbols are similar to those in Fig. 2.

8 Applications to a selection of solid rocky and ice-rich exoplanets

We applied the scaling laws derived here to a selection of exoplanets with masses comprised between 0.1 and 10 M_⊕, either a rocky or ice-rich composition, and a low surface temperature in order to avoid surface melting (for ice or rock; see Table 5). In order to derive a meaningful planet composition, we also chose planets for which both masses and radii were determined with relatively good accuracy (<5–10%). For each planet, we determined the composition in terms of iron content (for water-free rocky planets) or ice content (for ice-rich planets) compatible with the range of their observed mass and radius using the M–R relationship provided in Table 1. Subsequently, the Love number k₂ was determined using the scaling laws provided in Table 4. The Q factor was estimated from the values computed for planet masses of 1, 5 and 10 M_⊕ (see Figs. 7 and 10) and then extrapolated to the corresponding planet mass. As shown in Figs. 7 and 10, the Q factor strongly depends on the assumed viscosity. In Table 5, we provide the values estimated for a given viscosity profile in order to highlight the sensitivity to planet composition. For rocky planets, the Q factor is calculated by considering the reference viscosity profile (see Sect. 2) using Eq. (3). For ice-rich planets, the Q factor is estimated for a viscosity of high-pressure ice of 10¹⁶ Pa s (e.g. Durham et al. 1997), which can be considered as a reasonable lower limit (η_min). Table 5 and Figs. 12 and 13 summarize the results obtained for six exoplanets: three rocky ones and three ice-rich ones.

Fig. 10 Computed Q factor for ice-rich planets with masses of 1 (a), 5 (b), and 10 (c) M⊕, and three icecontents (10, 20, and 50%), as a function of viscosity η and tidal frequency χ for theice mantle. The vertical thick black segments indicate the expected values of Q for three planet candidates (Trappist-1 f, Trappist-1 g, and K2-240 c) for a high-pressure ice viscosity of 1016 Pa s. The black arrows indicate the trend if larger ice viscosity values are considered.

Fig. 11 Evolution of Q as a function of internal viscosity and tidal frequency for ice-rich and rocky planets.

8.1 Rocky planets

Kepler-36 b

This planet is consistent with an iron content larger than that of Earth, possibly up to 50%. However, due to the uncertainties on the mass determination, an iron content slightly lower than that of Earth (down to −23.5%) is also possible. Also, a small fraction of water cannot be excluded, but for simplicity we consider only a rocky composition.

LHS 1140 b

This planet is consistent with an iron content smaller than that of Earth possibly down to −50%. Here again, due to the relatively large uncertainty on the mass determination, a larger iron content cannot be excluded. Alternatively to a low iron content, the mass and radius of this planet could be explained by a small fraction of water ice (up to 5–10% depending on the iron content). We consider here only the rocky case.

Trappist-1 e

Among the three planets, this is the planet containing the largest iron content. Even if both radius and mass are determined with very good accuracy, the uncertainty on iron content remains relatively large as the M–R relationship is less sensitive to iron content in this mass range. As a consequence, a slightly subterrestrial iron content cannot be excluded. On the other hand, it is possible that the iron content exceeds 50% of the value for Earth. However, as we did not perform any computation above this value, we cannot estimate the upper limit with confidence in terms of iron content. This translates to an indetermination of the upper limit for k₂, indicated bya dotted line on Fig. 12.

The Love number, k₂, is mostly sensitive to the iron content of the planet, and also is slightly sensitive to planet mass for masses below ~2 M_⊕ (Fig. 12). Relatively large uncertainties on the iron content translate to equally large uncertainty on k₂. As already illustrated in Fig. 6, the Q factor is much less sensitive to iron content and depends mostly on the forcing period and the assumed internal viscosity. Asa consequence, assuming the same reference viscosity structure, the planet LHS 1140 b has the lowest Q factor among the three planets, as its orbital period is the largest one. As already mentioned, the values of Q factor listed in Table 5 have been calculated assuming the reference viscosity structure and an α parameter of 0.25. The viscosity may vary by a factor of 10–100 relative to the reference value, which results in Q values possibly three to four times larger. The values listed in Table 5 should be considered only as indicative. For any other viscosity and α values, the Q factor can be estimated using Eq. (3).

Fig. 12 Predicted k2 using Eq. (2) for a selection of three rocky exoplanets (Trappist-1e, Kepler-36 b, LHS 1140 b) from their mass, radii, and iron content, δFe derived from the M–R relationship provided in Table 1. Curve symbols are similar to those in Fig. 2.

8.2 Ice-rich planets

K2-240 c

We chose this planet as a possible representative of ice-rich planets of masses of the order of five times that of Earth. Among the catalog of detected exoplanets, we did not find any planet in this range of mass with surface temperature low enough to be compatible with cold icy surfaces. The surface temperature of K2-240 c is estimated to 389 K (Díez Alonso et al. 2018), which is the lowest temperature we found for this range of planet size. For this planet, only the radius has been constrained from transit measurements. The mass has been estimated using the only M–R relationship, following the law published by Weiss & Marcy (2014). Radial-velocity follow-up should be able to provide some constraints on the mass in the near future. Based on the published mass estimate, we concluded that this planet is an ice-rich planet with ice content approaching 50%, and possibly even more.

Trappist-1 f

This planet has an estimated surface temperature of 219 K, compatible with a cold icy surface. The mass and radius constraints imply that the planet contains at least 5% of water ice (assuming that the rocky part has a composition comparable to the Earth). Its ice content may reach values of 16%, and even larger values if iron-enrichment is considered in the rocky interior.

Trappist-1 g

This planet has an even lower surface temperature (199 K, Gillon et al. 2017) and ice content slightly higher than Trappist-1 f. Using our M–R relationship, we estimated the ice content to range between 13 and 26%.

The computed Love number for ice-rich planets is about half the values obtained for rocky planets at comparable masses (Fig. 13). By contrast, ice-rich planets are significantly more dissipative than rocky ones, with Q values possibly as low as 3; such as in the case of K2-240 c for example. This strongly depends on the assumed values for the viscosity of ice at high pressure however, which is not well constrained. The Q values listed in Table 5 correspond to the expected values for a viscosity of 10¹⁶ Pa s, which should be considered as the lowest estimate. For viscosity values ten times higher, Q can reach values exceeding 100 for Trappist-1 f and g, and thus become comparable to the Q factor predicted for rocky exoplanets. The Q factor provided here is only indicative and should be considered as the lowest estimate.

Fig. 13 Predicted k2 using Eq. (2) for a selection of three ice-rich exoplanets (Trappist-1 f, Trappist-1 g, K2-240 c) from their mass, radii, and the ice content, xice derived from the M–R relationship provided in Table 1. Curve symbols are similar to those in Fig. 2.

9 Conclusion

In the present study, we computed the tidal response of rocky and ice-rich solid exoplanets for masses ranging between 0.1 and 10 Earth masses. Accurate determination of the tidal parameters of exoplanets (k₂ and Q) taking into account their interior structure and rheology is essential to better characterize their rotational and orbital history, which has a direct impact on their climate stability and habitability. We showed that the amplitude of tidal response, characterized by the gravitational Love number, k₂, depends mostly on planet mass and the internal stratification (mostly controlled by the bulk composition). We found that for a given planet composition k₂ is mostly controlled by self-gravitation and increases as a function of planet mass; for a given mass, it depends on the relative size of the iron core, and hence on the bulk iron fraction. For a given mass, the presence of outer ice layers reduces the amplitude of tidal response compared to rocky planets with no ice. This reduction is controlled by the ice fraction. For both types of planet (rocky and ice-rich), we proposed relatively simple scaling laws to predict k₂ as a function of radius, planet mass, and composition.

For the Q factor, we did not find any direct control by the planet mass, which contrasts with the prediction of Efroimsky (2012) who assumes a homogeneous interior. However, as in this latter author, we observed a decrease of dissipation rate with increasing mass (larger Q with large planets), but not for the same reason. The reduction of dissipation rate is attributed here to the increase of mantle viscosity, which is expected to increase with pressure and hence with planet size (Tackley et al. 2013), while Efroimsky (2012) attributed it to self-gravity effects. While prediction using homogeneous body formulation is a relatively good approximation for small undifferentiated bodies, it does not describe correctly the tidal response of large planets. The effect of density stratification, the increase of elastic parameters and viscosity with increasing pressure, and the presence of a central liquid metallic core cannot be taken into account with such a simplified formulation.

Consistent with the results of Henning & Hurford (2014), we obtained a strong increase of dissipation rate (decrease of Q factor) for ice-rich planets compared to ice-free planets. This is explained by the relatively low viscosity of ice, which is expected to be several orders of magnitude below that of rocks. However, we should keep in mind that the effect of pressure is not well constrained and that it might approach values similar to rock near the melting point (10¹⁸–10¹⁹ Pa s). Viscosity is known to be strongly dependent on temperature, and therefore both temperature and viscosity profile should be modeled in a consistent way. In the present work, we assumed constant viscosity in each internal layer for simplicity and did not take into account coupling with thermal evolution as was done for instance by Běhounková et al. (2011) who dedicated their study to Earth-like planets. Future works are required to determine the influence of thermal structure on the dissipation rate and the retroaction of tidal dissipation on thermal evolution.

For planets orbiting at close distance to their stars, tidal dissipation may significantly contribute to the planetary heat budget. The amount of heat potentially dissipated depends primarily on the orbital configuration but also on the internal properties of the planets. The approach we developed in the present study allow the global tidal dissipation and its distribution in the planetary interior to be computed. We decided to focus here on the global tidal parameters, k₂ and Q, which are directly relevant to determining the retroaction of tidal friction on the orbital and rotational evolution of exoplanets. The distribution of heat dissipation in the interior and its implications for thermal evolution is addressed in a companion study dedicated to the Trappist-1 system (Breton et al. 2018, and in prep.). More generally, the methodology described in the present paper can be applied to any kind of solid exoplanet and can easily be coupled to any thermal or orbital evolution code.

Acknowledgements

The research leading to these results has received funding from the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007-2013 Grant Agreement no. 259285).

References

All Tables

All Figures

Fig. 1 Computed interior structure for a given mass and bulk composition. The bulk composition is defined relative to a reference composition corresponding to the solar value for the Fe/Si and Mg/Si ratio and to the terrestrial value for the Mg/(Mg+Fe) ratio in the silicate mantle. For rocky planets, the water ice mass ratio is 0%, while it can reach values up to 50% for ice-rich planets. For the examples shown here, an ice mass ratio of 20% is considered and planet masses of 1 and 10 M⊕ are displayed. The interior structure is divided into three main layers: ice mantle, rock mantle and iron core. The dashed lines indicate the separation between the low-pressure and high-pressure phases for both ice and silicate mantle.

In the text

Fig. 2 Mass–radius relationship, normalized to the Earth, obtained for rocky and ice-rich planets with various iron contents (−50% ≤ δFe ≤ + 50%) and ice fractions (0 ≤ xice ≤ 50%), relative to the composition of Earth. Each point corresponds to the results obtained using the modeling approach described above, based on Sotin et al. (2007). The lines correspond to the scaling laws using the polynomial coefficients provided in Table 1.

In the text

	Fig. 3 Comparison between our synthetic profiles (solid line) of density, ρ, bulk modulus, K, and shear modulus, μ, derived following the method described in Sect. 2 with those of the Preliminary Reference Earth Model (PREM; Dziewonski & Anderson 1981) (dotted line), assuming an interior composition corresponding to the reference case (see Table 2).
In the text

Fig. 4 Computed Q factor of the solid mantle of Earth for the 12.42-h tidal period using the synthetic profile displayed in Fig. 3 for various values of viscosity for the lower high-pressure (HP) mantle, ηHP (with ηLP∕ηHP = 100) and of α parameter for the Andrade rheological model. The bold black line indicates the estimated value of Q for the interior of Earth and the gray area shows the uncertainty range estimated by Ray et al. (2001).

In the text

	Fig. 5 Computed k2 for rocky planets with mass ranging between 0.1 and 10 Earth’s mass and iron content varying between −50 and +50% relative to the reference value, for a tidal period of 1 day. The lines represent the fit based on Eq. (2) and the polynomial coefficient provided in Table 4.
In the text

Fig. 6 Computed Q factor for rocky planets with mass ranging between 0.1 and 10 Earth masses (a) for a tidal period of 1 day and three iron content values varying between −50 and +50% relative to the reference value, and (b) for the iron content of the Earth and three tidal periods of 1 (square), 10 (triangle) and 100 (diamond) days. Two extreme values are considered for the very-high-pressure rock viscosity (ηVHP): 5 × 1022 Pa s (gray) assuming no viscosity increase relative to the reference Earth lower mantle and 5.1025 Pa.s (black) assuming a factor 1000 increase relative to the reference, i.e. the lower mantle of the Earth. Curve symbols are the same as in Fig. 2.

In the text

	Fig. 7 Scaling of the Q factor as a function of mantle viscosity η, tidal frequency χ, Andrade α parameter, and iron content, δFe, relative to reference values (ηref = 5 × 1022 Pa s, χref = 2π∕1day), for four planet masses (1, 2, 5 and 10 M⊕).
In the text

	Fig. 8 Computed k2 for ice-rich planets with masses ranging between 0.1 and 10 Earth masses and H2 0 mass fraction up to 50% relative to total mass. The lines represent the fit based on Eq. (2) and the polynomial coefficient provided in Table 4.
In the text

	Fig. 9 Computed Q factor for ice-rich planets with masses ranging between 0.1 and 10 Earth masses for two tidal periods (1 day in black and 100 days in gray) and three values of ice fraction (10% – cross, 20% – triangle, and 50% – diamond) for a HP ice viscosity of 1018 Pa s and thereference viscosity profile for the rocky part (see Sect. 2). Curve symbols are similar to those in Fig. 2.
In the text

Fig. 10 Computed Q factor for ice-rich planets with masses of 1 (a), 5 (b), and 10 (c) M⊕, and three icecontents (10, 20, and 50%), as a function of viscosity η and tidal frequency χ for theice mantle. The vertical thick black segments indicate the expected values of Q for three planet candidates (Trappist-1 f, Trappist-1 g, and K2-240 c) for a high-pressure ice viscosity of 1016 Pa s. The black arrows indicate the trend if larger ice viscosity values are considered.

In the text

	Fig. 11 Evolution of Q as a function of internal viscosity and tidal frequency for ice-rich and rocky planets.
In the text

	Fig. 12 Predicted k2 using Eq. (2) for a selection of three rocky exoplanets (Trappist-1e, Kepler-36 b, LHS 1140 b) from their mass, radii, and iron content, δFe derived from the M–R relationship provided in Table 1. Curve symbols are similar to those in Fig. 2.
In the text

	Fig. 13 Predicted k2 using Eq. (2) for a selection of three ice-rich exoplanets (Trappist-1 f, Trappist-1 g, K2-240 c) from their mass, radii, and the ice content, xice derived from the M–R relationship provided in Table 1. Curve symbols are similar to those in Fig. 2.
In the text