TIDAL HEATING IN MULTILAYERED TERRESTRIAL EXOPLANETS

Results are compared with a homogeneous tidal heating approach for exoplanets, described in detail in Henning et al. (2009). The problem of computing the tidal heat production rate for a simplified planetary body is approached using the viscoelastic form of the standard homogeneous tidal equation (Segatz et al. 1988):

$$\begin{equation} W_{{\rm tidal}} = -{\rm Im}(k_2)\frac{21}{2} \frac{R_{{\rm sec}}^5 n^5 e^2}{G}, \end{equation} \tag{ 1 }$$

where each parameter is given in Table 1, and Im(k₂) is the imaginary component of the complex valued load Love number. This term contains all information in this equation regarding viscous dissipation in the material, and is derived from a material's constitutive equation and complex valued rigidity. In Henning et al. (2009) Im(k₂) was computed for the Maxwell, Standard Anelastic Solid, and Burgers rheologies.

Table 1. Symbols and Parameters

Parameter	Symbol
Tidal heat rate	W
Semimajor axis	a
Eccentricity	e
Mass of the primary	Mpri
Mass of the secondary	Msec
Radius of the secondary	Rsec
Second-order gravitational Love number	k2
Second-order radial displacement Love number	h2
Second-order tangential displacement Love–Shida number	l2
Orbital mean motion	n
Orbital period	P
Pressure	p
Temperature	T
Ice layer thickness	tice
Displacement	ur, uθ, uϕ
Strain
Stress	σ
Internal gravitational potential	φ
External gravitational potential	Φ
Gravitational potential stress	ψ
Colatitude	θ
Longitude	ϕ
Colatitude relative to tidal symmetry axis	θT
Longitude relative to tidal symmetry axis	ϕT
Harmonic degree	ℓ
Density	ρ
Propagator matrix	Y
Aggregate propagator matrix	B
Propagator boundary matrix	M
Propagator solution vector	y
Propagator solution kernel	c
Propagator boundary vector	b
Seismic parameter	β
Quality factor	Q
Circularization timescale	τcirc
Maxwell timescale	τMax
Universal gas constant	R
Gravitational constant	G
Defining viscosity	ηo
Viscosity setpoint	ηset
Shear modulus	μ
Activation energy	E*
Activation volume	V*

Download table as:  ASCIITypeset image

The classical tidal circularization timescale may be written in terms of the energy loss rate to tides (Goldreich & Soter 1966; Murray & Dermott 2005),

$$\begin{equation} \tau _{{\rm circ}}=\frac{G M_{{\rm pri}} M_{{\rm sec}} e^2}{(1-e^2) a W_{{\rm tidal}} } \end{equation} \tag{ 2 }$$

where the (1 − e²) term can be neglected at small e. This form of the equation assumes negligible eccentricity change due to the tidal bulge raised on the primary by the secondary, which is true for most cases of a terrestrial class planet and star. This form of the circularization equation, in terms of W_tidal instead of in terms of Q, highlights the fact that circularization rates depend on a quantity that is time-dependent, and will vary both with the evolution of the system frequency n, and the internal planetary temperature.

Other forms of this equation expand W_tidal as in Equation (1), allowing cancelation of terms; however, this neglects verifying if a given planet is capable of supporting a particular dissipation rate. If during circularization, large-scale melting occurs, then W_tidal will change. In many cases, melting will lead to piped advection of melt (O'Reilly & Davies 1981; Moore 2001; Monnereau & Dubuffet 2002; Connolly et al. 2009; Spiegelman 1993; Ribe 1987) either to a surface ocean, subsurface ocean, or both. Addition of an ocean to an otherwise initially solid world can have a significant impact of the total tidal dissipation, and thus alter τ_circ, as discussed in Section 5.

Quality factors, introduced from seismology, are a useful although occasionally misleading aggregate scalar factor to describe inverse energy loss per cycle of a damped oscillator. Q factors are often treated as effectively constant, although in reality they are strong functions of forcing frequency, temperature, layer structure, and numerous internal material properties that vary with time. Values of Q are related to the lag angle of a tidal bulge _t via the formula: Q = sin |_t|. For Earth, a Q of 12–34 has been measured based on the expansion of the moon's orbit (Yoder 1995; Dickey et al. 1994; Murray & Dermott 2005); however, such low values (thus high dissipation) are due to the lunar-tide induced flow in Earth's water ocean. Goldreich & Soter (1966) estimated Q for Mercury, Venus, and the Moon as Q_merc ⩽ 190, Q_venus ⩽ 17, and 10 ⩽ Q_moon ⩽ 150. However, these values for Venus and Mercury make several assumptions about each world's spin history. Lainey et al. (2007) recently estimated Q_mars = 79.91 (±0.69) from the motion of Phobos. Data summarized in Karato & Spetzler (1990) suggest Earth's lower mantle Q at tidal periods from 1–10 days is in the range 50–200. Based on such results, a common assumption for a dry 1 M_E exoplanet has been the round number 100. In reality, Q is coupled in a viscoelastic system with the second order tidal Love number k₂ (Segatz et al. 1988), leading to the value Im(k₂) in Equation (1), which replaces the ratio k₂/Q in the non-viscoelastic form of Equation (1) (e.g., Peale et al. 1979). For use in common astrophysical equations, however, it is often useful to extract from Im(k₂) an effective Q value, by the formula Q_eff = Re(k₂)/− Im(k₂). Earth-sized effective Q values range as far as ∼1–10⁷.

For solid bodies, the common practice of adopting a fixed Q factor and calculating a resulting circularization time can lead to completely incorrect predictions, because it skips the step of examining if the tidal heating rate of the planet is realistic or sustainable. Consider GJ 876 d, with a minimum mass of 6.8 M_E, period of 1.937 days, and host mass of 0.33 M_sol. Using a fixed Q of 100 for this world leads to a rapid circularization timescale of 4 million years. This, however, would require the planet output 80 million TW of power for all 4 million years, which is extremely unlikely (the global heat rate of the modern Earth is only ∼44 TW). Instead, the planet will experience global-scale partial melting, altering the initial assumed Q. Using a tidal equilibrium model (Section 3.3) with partial melting, and a starting eccentricity of 0.139 (Correia et al. 2010), leads to a valid sustainable tidal dissipation rate of only 80,000 TW, equivalent to a global effective Q value of 100,000, and a circularization timescale of 4 billion years. A main argument to proceed with the research in this work is the primary concern of fixed Q models for icy or rocky planets, namely that Q is highly sensitive to internal structure, and the need to assess how Q (or by proxy Im(k₂)) varies when global-scale partial melting or liquid magma or water ocean layers exist. This work may be seen as a precursor to the larger goal of assessing the tidal evolution of systems including terrestrial planets in time, particularly during orbital scattering and circularization, when such ocean layers (or simply hotter and less viscous layers) are very likely to develop.

Our goal is to understand the grand scale differences in the response of planets of different bulk structures, at a time when even the rheological details of Earth's own deep interior remain in debate. For bulk tidal dissipation as a function of depth, the critical parameters are the viscosity, shear modulus, and density, with viscosity by far exceeding all other parameters in terms of control of the final tidal outcome. This dominance of viscosity allows one to begin with several simplifying assumptions in regards to density, rigidity, and composition.

For a homogenous Maxwell rheology planet, −Im(k₂) in Equation (1) may be expressed analytically (Henning et al. 2009)

$$\begin{equation} -{\rm Im}(k_2)_{{\rm Max}} = \displaystyle\frac{57 \eta \omega }{4 \rho g R_{{\rm sec}} (1+ [ (1 + (9.5 \mu \rho ^{-1} g^{-1} R_{{\rm sec}}^{-1}))^2 \eta ^2 \omega ^2 \mu ^{-2} ])}, \end{equation} \tag{ 3 }$$

where η is viscosity (in Pascal-seconds), μ is the shear modulus (sometimes referred to as rigidity), ω is the forcing frequency, ρ is the homogeneous planet density, g is surface gravity, and R_sec the planet radius. Similar expressions may be written for other rheologies (Henning et al. 2009). Note that Equation (3) embeds an assumption of incompressibility by considering only deviatoric stress components within a world, and can thus become less accurate for large worlds such as mini-Neptune class planets.

While it is possible to write the analytical expression of Im(k₂) for multilayer planets, even for two-layer bodies such expressions rapidly become impractical to work with. For several layers, they readily require several pages to express. We therefore shift to staged numerical computation to determine Im(k₂) in multilayer cases. The essential character of Equation (3), however, remains the same, in that outcomes are a result of η, μ, and ρ for each layer.

For the shear modulus of silicates we adopt a baseline value of 5 × 10¹⁰ Pa unless otherwise noted, such as in some Earth models below. For solid ice, we adopt a standard value of 4 × 10⁹ Pa. While other works use a variety of shear modulus values for both material classes, what is most important (especially considering the unknown role of impurities) is the comparatively very small range of uncertainty in this parameter relative to uncertainty in viscosity. In general, impurities weaken a pure crystalline material, and therefore the inclusion in water ice of other volatiles such as CO₂, NH₃, and CH₄, either as clathrates or as ices themselves, may be considered in a simplified fashion by considering lowered μ values. Because of high pressures, weakening due to faults and material damage is not a major consideration below the first 1–2 km of a major planet.

Viscoelastic materials typically experience a material resonance as a function of forcing frequency (see e.g., Nowick & Berry 1972), at which maximum dissipation occurs. For the Maxwell rheology, this peak response occurs when a forcing period equals the Maxwell timescale τ_Max = η/μ. A multilayer planet will not have one characteristic Maxwell timescale but many, and different layers will resonate at different forcing rates. Typical Maxwell times for silicates are shown in Figure 1 and cover a broad range from hours to years. For typical ice with a viscosity of 10¹³, 10¹⁴, and 10¹⁵ Pa s, τ_Max = 0.6 hr, 6.9 hr, and 2.8 days. The relative separation between material layer Maxwell times and forcing periods plays a central role in determining tidal outcomes. Layers that are viscoelastically well-matched to their forcing periods generally experience the highest tidal heating rates. For silicates the very wide distribution of η and μ values makes estimation of behaviors difficult prior to more individualized analysis. For ices, most cooler material will have an optimal response in the lower typical range of planetary and exomoon periods.

Zoom In Zoom Out Reset image size

For silicates and ices, shear moduli are nearly constant with temperature up until very near the melting temperature, at which point various models of shear weakening begin to occur (Fischer & Spohn 1990; Berckhemer et al. 1982). For pure water ice with a single melting temperature, the shear modulus falls directly to zero upon melting. Silicates melt across a range of temperatures, beginning at the solidus T_sol, with zero melt fraction, and ending at the liquidus T_liq with a melt fraction of 1.0 (with melt fraction varying linearly in temperature in between). A typical surface range for silicate melting is T_sol = 1600 K, and T_liq = 2000 K (Henning et al. 2009). Higher pressure melting points may be found via the Simon equation (Poirier 2000). The shear modulus for silicates remains close to its solid value for approximately half the melting range, and falls rapidly to zero at the point where individual mineral grains loose physical contact with one another, within a surrounding liquid bath. This point is referred to as the breakdown (or disaggregation) temperature T_br, and typically lies between 40% and 60% melt fraction (Moore 2003b).

Density may be treated either by utilizing fixed average bulk values appropriate for a given layer, or through detailed equation of state calculations, which in turn depend on a temperature profile. Density primarily influences the elastic component of the load Love number (e.g., Re(k₂)) for a planet on which tidal heating depends only linearly. Such Love numbers for the objects studied here vary typically by less than a factor of ∼2 due to the large overall mass. By contrast, viscosity variations may readily alter heat outputs by 4–5 orders of magnitude. Therefore, for simplicity we focus on cases with layer densities assigned (or in the case of Earth, taken from measurements) rather than calculated from equations of state and accompanying temperature/convection assumptions.

The tidal response of a planetary iron core is often weak, however, it is straightforward to include the core in multilayer calculations for completeness, and to verify such weakness, especially since it comes at little computational cost. Data on the viscosity of core analog Fe or Fe Ni alloy materials is limited. Frost & Ashby (1982) use a value of 1 × 10¹³ Pa s. Dumberry & Bloxham (2002) use a range from 1 × 10¹¹ to 1 × 10²⁰ Pa s. Buffett (1997) suggests a value less than 1 × 10¹⁶ Pa s. We discuss the results of such uncertainty further in Section 5.

Liquid metal in an outer core has very low viscosity. Ab initio methods by de Wijs et al. (1998) suggest values of 1.3 × 10⁻² Pa s for an inner core boundary at 6000 K, 1.2 × 10⁻² Pa s for a core mantle boundary of 4300 K, and 1.5 × 10⁻² Pa s for a core mantle boundary of 3500 K. Earth observations from the Chandler wobble and free oscillations however lead to high values of 10⁴–10¹¹ Pa s (Verhoogen 1974; Won & Kuo 1973) for the outer core. However, since even the highest value in this range is two orders of magnitude below the lowest solid viscosity used, tidal heating in any outer core is always found to be entirely negligible.

Viscosity is strongly dependent upon numerous material parameters, the most important of which is temperature. Pressure, grain size, and the total stress magnitude also play central roles. Stress state is most important for determining which modes of viscous behavior dominate the strain rate for an applied loading. Microscale viscous behavior may occur via diffusion creep (where crystal lattice voids migrate within grains), dislocation creep (the migration of linear lattice nonconformities known as dislocations), slippage along grain boundaries, or the diffusion of grain boundaries. Of these mechanisms, Moore (2006) demonstrates how diffusion creep dominates for ice under a wide range of tidal conditions and grain sizes.

As with shear moduli, in this work we focus on selecting informed estimates for viscosity in each material layer, as opposed to deriving such viscosities directly from a temperature and pressure profile. There are several reasons for this choice. First, the method of deriving viscosities from first principles at depth in a planetary body has not yet been fully successful even for the Earth, where postglacial rebound studies have provided measurements of the full mantle viscosity profile (Mitrovica & Forte 2004), but matching such values requires numerous posteriori adjustments to theoretical models, such as invoking the lateral inhomogeneity of rising and falling convection plumes and uncertainties in petrological composition. In reality, is it likely that the pressure dependence of planetary materials as currently modeled for low pressures does not extrapolate well to deep mantle pressures. For the Earth, rising temperatures act to lower viscosities, while high pressures act to raise viscosities, and results are extremely sensitive to which of these two processes is allowed to dominate, in effect meaning results are largely controlled by the choice of the activation volume V* (see Equation (4)), a value for which there is little laboratory data at high pressures.

In addition, deriving viscosities from a temperature and pressure model by depth inherently converts the problem into a highly unstable iterative computation. Two strong feedbacks are intertwined in this case. The vigor of convection, which determines upper boundary layer thicknesses and heat flux rates, is itself a strong function of viscosity. Second, heat flux rates depend upon the level of tidal heating, and thus tides strongly control planetary temperatures and viscosities. In reality, these feedbacks are resolved in time for a given orbit by an equilibrium between tidal heat input and convective heat loss. Planets will often evolve in time (typically a few million years) toward a very stable equilibrium temperature (Moore 2003b) where tidal heat input balances convective heat loss. In a large fraction of rheological and orbital cases (Henning et al. 2009), terrestrial exoplanets will reach tidal-convective equilibrium very close to either the solidus T_sol or breakdown T_br temperatures. Therefore, for high tidal-heat Hot Earth cases, the most important viscosities to consider are those expected to reside near these states.

Previous analysis has mainly considered such tidal-convective equilibration for mantles characterized by a single bulk temperature and single melt fraction. While this is not a bad assumption for a well-mixed adiabatic convecting mantle, it may become significantly inappropriate once large-scale melting and melt mobilization begins for hotter planets. Multilayer models, where liquid or weak solid partial melt layers may be explicitly introduced along with homogeneous convecting layers, offer a significant improvement in their ability to determine a true range of tidal outcomes.

For silicates, viscosity prior to melting may be computed from an Arrhenius model

$$\begin{equation} \eta = \eta _o e^{\frac{E^{\ast } + pV^{\ast } }{RT} }, \end{equation} \tag{ 4 }$$

where R is the universal gas constant, E* is an activation energy, T is temperature, p is pressure, and V* is an activation volume. A defining viscosity η_o coefficient is calculated based on a viscous setpoint, for example, η(T = 1000 K) = η_set. In Figure 1 we use η_set = 1 × 10²⁰, 1 × 10²¹, 1 × 10²² Pa s for a weak, moderate, and strong silicate rheology, and compute a corresponding η_o for each. The low value for defining viscosity represents a more volatile-rich mantle (Hirth & Kohlstedt 1996), while the high value represents a more devolatilized case. Viscosity reductions due to temperature and melting occur on top of this setpoint selection. Therefore, η_set is a compositional choice, and is not the actual viscosity of a planet's mantle, which is likely to be much lower, due to a significant partial melt fraction.

A complex falloff of viscosity occurs for silicates near the solidus temperature (Moore 2001, 2003b; Fischer & Spohn 1990; Berckhemer et al. 1982). We omit the details of such models here, except to note that the minimum viscosity for silicate materials prior to complete melting is in the range 1 × 10¹²–1 × 10¹⁶ (a range notably similar to the viscosity of ice).

The majority of tidal dissipation on ice–silicate hybrid planets occurs within ice layers. For ice, viscosity prior to melting may be computed following Goldsby & Kohlstedt (2001),

$$\begin{equation} \dot{\epsilon } = A_x \sigma ^{n_{\sigma }} h^{-m_g} e^{ - \frac{E^{\ast } + pV^{\ast } }{RT} } \end{equation} \tag{ 5 }$$

where $\dot{\epsilon }$ is strain rate, A_x is a creep mechanism scale parameter, σ is stress, and h grain size. For volume diffusion creep, the exponent values are set at n_σ = 1 and m_g = 2, and we may follow Moore (2006) for Ice Ih to adopt A_x = 9.06 × 10⁻⁸/T Pa$^{-n_{\sigma }}$ m$^m_g$ s⁻¹, E* = 59,400 J mol⁻¹, and V* = −1.3 × 10⁻⁵ m³ mol⁻¹. Typical values for h include the range 0.0001 to 0.01 m. Note that the activation volume here is negative, and therefore Ice Ih is expected to decrease in viscosity at higher pressures, unlike silicates that rise in viscosity. Figure 2 shows the results of Equation (5) for a range of pressures and temperatures appropriate for an exoplanet ice mantle. Typical pressures for such ice mantles are demonstrated in Figure 3. The pressure variation must be considered cautiously however, because the activation volume for high pressure ice phases, especially Ice X, are not well known. Figure 2 does, however, allow us our goal of a range of informed estimates for ice viscosities for the cases of interest, even if not dynamically linked by all feedback loops to the real temperature-depth profile of each world.

Zoom In Zoom Out Reset image size

The compositional complexity of ice can be great, both from the phase diagram of water ice, as well as from the inclusion of multiple ice species such as CH₄, NH₃, and CO₂. Eutectic blends of water with ammonia are often invoked as a likely mechanism for melting point depression in outer icy satellites. As with shear moduli, it is better to encompass such complexity by testing a range of candidate ice viscosities. This testing method has the additional benefit of reducing cases, since a very large range of temperature and compositional realities may eventually lead to the same viscosity, a handful of tests over the full viscosity range may instead then be reverse-mapped back to a heavily degenerate range of possible material causes.

A common central value for ice viscosity for outer moons is 10¹⁴ Pa s (Poirier et al. 1981; Goldsby & Kohlstedt 2001). The temperature-only-dependent range of Equation (5) in Figure 2 is 1 × 10¹⁷ to 1 × 10¹³ Pa s, and we therefore tested a range of values at least two orders of magnitude around a baseline of 1 × 10¹⁵ Pa s, with low values best representing compositional impurity, and higher values better representing cold ices.

Before extending an Earth model toward a super-Earth model, we first seek to establish and characterize a reliable multilayer Earth model. In previous work we have studied homogeneous tidal models of extrasolar Earth mass planets. Thus it is also of central value to investigate the degree to which a range of improved multilayer Earth models alters the outcome of tidal heating.

We find that an essential feature required to reproduce an Earth-like Re(k₂) value of around 0.299, is a high value for the shear modulus of the silicate mantle. A typical value for the rigidity of surface rocks is 5 × 10¹⁰ Pa. For the deep interior, seismic velocity studies provide the parameter β = (μ/ρ)^1/2, which may be inverted to determine estimates of μ given an estimate for the density. Lower mantle densities are expected to vary in the range from 5560 to 4380 kg m⁻³. A nominal value of 4600 kg m⁻³, coupled with an outer core density of 10,200 kg m⁻³ and inner core density of 12,000 kg m⁻³ leads to a nominal ∼1 M_E total mass, and to shear moduli in the range from 1.6 × 10¹¹ to 2.4 × 10¹¹ Pa. An adopted value of μ = 2.25 × 10¹¹ Pa reproduces the expected Re(k₂) = 0.299 for such a three-layer model.

For a four-layer model, we divide the mantle into an upper and lower layer at the 670 km depth mark. Using an upper mantle estimated density of 3900 kg m⁻³, a lower mantle density of 5000 kg m⁻³, upper mantle rigidity of ∼7 × 10¹⁰ Pa (from β of the upper mantle), increases the lower mantle rigidity required to achieve the expected Re(k₂) up to 2.55 × 10¹¹ Pa. Regardless of these specific values, the primary point is that using a surface rigidity of ∼5 × 10¹⁰ at depth can lead to unrealistically high Re(k₂) values of ∼0.6–0.7.

To model the continuous variation of properties within the Earth we implement a more detailed model using the Preliminary Reference Earth Model (PREM; Dziewonski & Anderson 1981). This model includes 33 layers, at roughly 200 km steps, and resolves Earth's crust, Low Velocity Zone, Upper Mantle Transition Zone, and D'' layers. For computational reasons the liquid outer core is still treated as one layer with an average bulk density. Following PREM, β in the lower mantle lies in the range 5945.1 to 7264.7 ms⁻¹. Mapping these β values against the range of densities leads to shear moduli from 1.5 × 10¹¹ Pa to 2.9 × 10¹¹ Pa, all much greater than the surface value of 5 × 10¹⁰. Repeating this exercise for the upper mantle yields β = 5570.2–4491.0 ms⁻¹. Because viscosities are not resolved by the seismic input data of PREM, in Figure 4 we test both a version of PREM with constant viscosities throughout the mantle, as well as a model with viscosities adapted from post-glacial rebound studies (Mitrovica & Forte 2004).

Zoom In Zoom Out Reset image size

For the model and forcing range of primary interest, the selection of viscosity has a weak interaction with Re(k₂), but is the principal control of energy dissipation via Im(k₂). In Equation (1), −Im(k₂) is effectively used as a replacement for Re(k₂)/Q. Thus we seek viscosity values which bring a baseline Earth model's −Im(k₂) close to that required to achieve a reasonable global effective Q_eff value. This, however, is complicated by the strong frequency dependence of Im(k₂), and on the fact that we are modeling the solid tidal dissipation of worlds, while the Q of the Earth from lunar retreat rates, Q ∼12, is mainly controlled by wave activity in Earth's surface water oceans. Even when we include water oceans in models in this work, we do so for their decoupling impact on layers, and we are not including wave dissipation or dissipation due to flow in or around specific geometries (Tyler 2008; Chen et al. 2014). The quasi-static contribution of water oceans to dissipation is otherwise negligible. The Q value for a dry Earth is unknown. There are two forms of dry in this context: ocean-free, and free of large-scale hydration in mantle minerals. A planet with little or no mantle mineral hydration may be highly viscous, as H₂O can lower viscosities by several orders of magnitude (O'Connell & Budiansky 1977; Hirth & Kohlstedt 1996). For ice–silicate super-Earths, however, this is unlikely, given that a large ice mantle implies a high water fraction to begin with, thus we are interested in a Q for a baseline iron-rock core to these hybrid worlds that is free of the dissipation contribution of a surface water ocean, but is not altogether dry mineralogically. Q = 100 is often used in the literature as a dry 1 M_E planet value, and the authors have used Q = 50 previously as a compromise for ocean free but still hydrated 1 M_E worlds.

Bulk viscosities around 1 × 10²⁰ Pa s for the rigidity model previously cited, at a period of 30 days leads to Q_eff = 1445, and at P = 300 days to Q_eff = 145. Values for silicate mantle viscosity around 1 × 10²¹ Pa s typically push Q values above 10,000 as in Figure 5. At a 30 day period, bulk viscosities of 1 × 10¹⁹ Pa s can lead to Q_eff = 149, 5 × 10¹⁸ Pa s to Q = 75, and 1 × 10¹⁸ to Q = 15. As seen in Section 4.2, as soon as an ice mantle is applied to a world, it is the ice viscosity and not the silicate viscosity that dominates the planetary tidal result.

In Figure 4 we model multiple variations in Earth structure. In each case the total dissipation is reported in the legend, and is a function of the applied e = 0.1, P = 30 days, and 1 M_Sol orbital model. Case 4A is a homogeneous model, set to the average bulk density of the Earth, with a shear modulus of 5 × 10¹⁰ Pa and viscosity of 1 × 10²¹ Pa s, following previously published baseline homogeneous Earth analogs in Henning et al. (2009). In homogeneous cases generally, dissipation is maximum in the mid-mantle, finite at the surface, and zero at the core.

Zoom In Zoom Out Reset image size

Note that text references to lettered Cases in figures are prefixed by figure number to help reduce confusion due to the large number of possible cases to be shown.

In Case 4B, only a solid core and solid mantle exist, to model a simple two layer world, equivalent to a highly evolved Earth where the full core has crystallized. A liquid outer core will decouple the inner core from major tidal flexure, such that inner core dissipation is generally weak, and the viscosity selected for the inner core is of low importance. When the solid silicate and iron layers are directly coupled, the iron viscosity is of central importance. Low values around 1 × 10¹⁷ Pa s used elsewhere lead to core dissipation dominating the response by a full order of magnitude. However, since the viscosity of such an iron core is highly uncertain, especially due to the uncertain geotherm of such a high age world, we instead plot in Figure 4 a case where the viscosity of both the iron and mantle are matched at 1 × 10²¹ Pa s each, in order to detect the difference due to the density and shear moduli otherwise hidden by a viscosity contrast.

Case 4C includes a liquid outer core where dissipation is negligible. This decoupling allows greater flexure of the mantle and thus greater dissipation. Inner core viscosity is switched to 1 × 10¹⁷ Pa s in this case, while mantle viscosity remains at 1 × 10²¹ Pa s. Case 4D begins to approach a realistic Earth model, with an inner and outer core, upper and lower mantle, and low dissipation lithosphere.

The final two Cases 4E and 4F invoke the PREM model through 33 material layers where density and shear modulus vary continuously. The PREM model includes the D'' layer at the core-mantle boundary (CMB), which is found to exhibit high dissipation due to its own low strength blending with the natural maxima of dissipation at the CMB from Cases 4C and 4D. Extra layers were added at this location to enhance resolution. Tidal enhancement in the upper mantle and attenuation in the lithosphere for each PREM model is similar in form to Case 4D where properties had remained constant. The PREM model does not describe viscosities with depth, therefore Cases 4E and 4F test two viscosity models. In Case 4E, viscosities are roughly based on a synopsis of models from Figure 1 of Mitrovica & Forte (2004), where high viscosity occurs both at high depth due to high pressure, and at shallower depths due to lower temperatures, with lower values in the D'' region, mid-mantle, and upper mantle. The modeling here is not intended to be exact, but to illustrate the impact of such viscosity variations on tidal output. In Case 4F, constant viscosities of 1 × 10²¹, 1 × 10²⁰, and 1 × 10²¹ are applied through the lower mantle, upper mantle, and lithosphere, such that the impact of the PREM gradients in density and shear modulus relative to Case 4D can be seen.

Figure 4 does not include a model with a very weak asthenosphere or magma ocean (Tonks & Melosh 1993; Elkins-Tanton et al. 2003; Elkins-Tanton 2008), however, such models were also tested. Inclusion of a magma ocean just below the lithosphere often has negligible impact on the dissipation of lower layers, and the high rigidity of the lithosphere typically continues to generate negligible heating despite this decoupling. The primary impact of magma ocean inclusion is therefore the simple removal of upper mantle dissipation throughout the volume assumed to instead be liquid. Inertial effects, waves, and tidal flow (Tyler 2008) in a magma ocean, not modeled here, may however contribute in a similar manner to Earth's water ocean.

The primary result from Figure 4 is that non-homogeneous cases often produce tidal dissipation greater than the homogenous model, and very rarely produce less. Further enhancement may occur with large weak asthenospheres as in Figure 5. Enhancements by a factor of ∼1.25–1.50 are common. This occurs for two main reasons: first, including a liquid outer core allows greater tidal flexure, and second, more detailed models typically allow inclusion of local low viscosity zones. These effects often dominate over the reduction in heating caused by the removal of the liquid outer core volume from contributing to global heating. Even when low viscosity zones are small compared to overall planet volume, their dissipation typically dominates, similar to the expected dominance of the asthenosphere of Io (Segatz et al. 1988). The relative shape of depth profiles for the above cases at P = 2 or 365 days are similar to those in Figure 4, despite dramatically higher/lower overall magnitudes. At P = 3 days, total dissipation for Cases 4A–4F are: 96, 185, 281, 295, 172, 247 TW, or enhancements by factors of ∼1.8–3.1, with all multilayer results higher than the homogeneous case. Values for alternate eccentricities or host masses may be scaled using Equation (1).

Zoom In Zoom Out Reset image size

Figure 5 shows the distribution of tidal heat for a multilayer Earth, both in depth and across the surface in latitude and longitude, with a focus on variations in asthenosphere viscosity. Moderately low viscosity silicate asthenospheres generally attain the highest tidal response. Upper mantle viscosities are selected here to show the transition between two common surface patterns. Densities for this model are 12,000, 10,000, 5000, 3900, 2600 kg m⁻³ for the inner core, outer core, lower mantle, upper mantle, and lithosphere. Shear moduli are: 1.5 × 10¹¹, 0, 2.5 × 10¹¹, 7 × 10¹⁰, 5 × 10¹⁰ Pa. The inner core viscosity for this model is 1 × 10¹⁶ Pa s and the lithosphere viscosity 1 × 10²¹ Pa s. The top left panel of Figure 5 is remarkable as it is the inverse of the asthenosphere dominated six-lobed pattern found for Io (Segatz et al. 1988). We find that the mid-valued viscosity contrast pattern (middle panels) is most common for an Earth analog exoworld.

Variation of the orbital period has minimal impact on the pattern of the surface distribution when outside of the transition range where upper and lower mantle heat rates begin to converge, but can dramatically change the pattern when this range is encountered. For most models, a period of 3 days leads to an upper mantle dominated pattern, while a 300 day period is less well tuned, causing a decrease in heat from lower viscosity zones, switching behavior to a lower mantle dominated pattern. Patterns of surface heat flux are often transitional for periods from 10 to 50 days.

We may roughly divide ice–silicate hybrid exoplanets into those which may contain internal liquid water oceans, and those which do not. Water oceans are expected to be a shallow or surface phenomenon (Fu et al. 2010), while in contrast, ice mantles may plausibly exist with very large thicknesses, up to the limit of the ∼2 R_E (∼12,700 km) ice mantle thickness expected for 3.8 R_E Neptune. In practice we test the impact of internal water oceans within ice mantles of up to 3000 km thickness, and assume that above this ice thickness, basal or internal oceans become unlikely and interest lies mainly in the physics of the solid ice layers fully coupled to solid silicates below. We first consider ice-only mantles, and then consider ice mantles in the presence of water oceans. In order to focus on bulk structural questions first, in this section we report upon ice mantles with constant viscosities and shear moduli, however, in Section 5 we consider gradations in material properties.

Figure 6 shows the distribution of tidal heating with depth for a range of ice–silicate hybrid planet models. Two key results are found: first, that ice–silicate tidal systems are highly insensitive to the choice of silicate/iron core model, and second, that the primary controls on the magnitude of total tidal heating (and thus effective Q values) are the viscosity of ice used, and the presence or absence of any decoupling water ocean.

Zoom In Zoom Out Reset image size

Iron-core mass fraction changes are found to have only a very minor impact on overall results due to the dominance of ice material properties on the tidal response. High iron core fractions are similar in impact to higher silicate densities in the sense that the entire iron-silicate component of the planet acts as an aggregate core. Switching from a liquid iron outer core, to an all-solid or all-liquid iron core has limited impact beyond planets with very thin ice mantles for the same reason.

Surface maps of ice-hybrid super-Earth tidal heat flux integrated over one orbit and summed over depth are shown in Figure 7. Ice mantles in general often produce response patterns in latitude and longitude similar to the asthenosphere-dominated response described in Segatz et al. (1988) for Io. In general, cases involving any low viscosity material of modest depth overlying a large depth of high viscosity material generate six-lobed global maps of dissipation with maxima at or near the equator. In contrast, when dissipation is optimal in a material which spans the majority of the body, as in the mantle-dominated cases of Segatz et al. (1988) for Io, or any homogeneous planet model, then tidal heating is maximum at the poles with two broad minima at the equator. Such polar maxima for tidal heat production may roughly be understood to arise from the dominance of shear stress terms in the overall system response.

Zoom In Zoom Out Reset image size

A new situation occurs in Figure 7, panel 5, where the thickness of the low viscosity material with the more optimal tidal tuning is no longer a thin top veneer, but instead approaches one-third of the body radius. In these situations the pattern of dissipation evolves beyond the six lobed case of a thin active upper layer, and heat becomes further concentrated to a belt around equator. When this layer approaches half of the body radius, as in Figure 7, panel 6, or the top panel of Figure 5, the pattern transforms to one of six minima, or the inverse of Figure 7 panel 3. Note that previous work has described the range in outcomes from Figure 7 panels 1 to 3 as endmembers of the solution space, however Figure 7 demonstrates that these patterns are in fact part of a much larger system of solutions.

The transition in surface patterns shown in Figures 5 and 7 is partly explained via Figure 8, which demonstrates how tidal heating as a function of depth systematically changes in an low-viscosity upper layer of increasing thickness. Homogeneous ice mantles of viscosity 1 × 10¹⁵ Pa s, density 1000 kg m⁻³, and shear modulus of 4 × 10⁹ Pa are shown from 400 km to 10,000 km thickness. In all cases the silicate/iron core is a homogenous solid with constant density set to match 1 M_E, because of the demonstration in Figure 6 that ice mantle tidal results are insensitive to silicate/iron core structure, especially for ice mantles above a few hundred kilometers thickness. Again, an orbit of e = 0.1, P = 60 days at an analog of 55 Cnc is used, thought the phenomenon highlighted here is not sensitive to these choices, only the total heat magnitudes are.

Zoom In Zoom Out Reset image size

As an ice mantle grows in thickness, it transitions from a profile where tidal heating is maximum at the base of the layer, to a form where a mid-layer maxima rises to dominate. Note that the mid-layer feature need not actually exceed the basal maxima in order to dominate the result, because the feature's widths also matter when total power is integrated with depth (Note that in units of W m⁻¹, the higher importance of upper layers due to their greater volume has already been taken into account). Results in Figure 7 may thus be reinterpreted as follows: at t_ice = 100 km, silicate-core features remain the dominant surface expression. At t_ice = 400 km, ice heating features dominate, with a basal maxima in layer heating leading to a six-lobed surface pattern. At t_ice = 3000 km, the six-lobed basal heating feature begins to be replaced by the mid-layer maxima, with the net result of an equatorial focusing of heat. At t_ice = 6000 km, the mid-layer maxima in heating fully dominates surface expression. Above this thickness the top layer is sufficiently thick that the surface expression begins to return to the form of a homogeneous body, with the effect of any small poorly coupled core eventually vanishing in importance. We find this pattern cycle is common for any high tidal susceptibility material placed atop a low tidal susceptibility core, and considered in terms of increasing thickness (e.g., it also applies to a growing silicate asthenosphere).

The origin of the multiple tidal heating surface patterns which may occur has no simple origin, and arises from the summation of dissipation maps which can be decomposed by tensor product component as in Figure 9. Such component maps may also be considered for the instantaneous work in time at each orbital position, as well as for the stress and strain tensor terms themselves. Individual maps vary with depth and forcing frequency, such that their collective sum weighed by varying magnitude in depth becomes more difficult to visualize. For two-layer worlds, Beuthe (2013) discusses the origin and range of tidal heating patterns in detail. In general, shear terms in stress lead to greater polar heating, and are dominant for the typical mid-layer maxima of homogeneous planets or extremely deep asthenosphere-like layers. The basal maxima typical of the CMB or any ice–silicate transition contain a more even mix of shear and normal stress terms, leading to more equatorial final expressions of surface heat.

Zoom In Zoom Out Reset image size

There is sufficient systematic variation in tidal surface patterns based on internal structures, that they may help to constrain the interiors of observed exoplanets in the future. Observations are already able to map large-scale hemispheric temperature dichotomies and offsets in temperature maxima from the substellar point for Hot Jupiters (Knutson et al. 2007; Cowan & Agol 2008). Studies suggest that mapping Earth analog solid surfaces due to rotation in inclined orbits may be plausible with upcoming technology, including both 1D maps (Fujii et al. 2011), and 2D maps (Kawahara & Fujii 2010, 2011; Fujii & Kawahara 2012). To observe the heat pattern of a super-Earth due to its internal tidal activity would require a planet orbiting close a very low luminosity host, ideally in a transiting configuration which allowed a tight constraint on orientation of the equator. Large moons far from host stars may be advantageous targets, and Peters & Turner (2013) have discussed the possibility of directly imaging the tidal signature of such hot exomoons. A rare ideal case would include periodic eclipse of one pole only, allowing for a polar vs. equatorial relative measurement. If the lightcurve of a tidal planet could be shown to have four brightness maxima in longitude, this would indicate a thin tidally active upper layer, but only two maxima would point toward a tidally active layer spanning the majority of the radius. Liquid oceans and atmospheres however will redistribute heat and may blur such effects. Yet even with the complexity evident in observing Io, or the lack of global distribution of heat on Enceladus, patterns of tidal enhancement have helped to constrain interior models (Kirchoff et al. 2011; Hamilton et al. 2013), and any opportunity to constrain the interior of a world (Ragozzine & Wolf 2009) in another starsystem will be significant.

Figure 8 also includes the total heat output and effective Q values for worlds with deep ice mantles. Such total heat values present a mixed picture. On one hand, values are mostly small compared to the ∼44 TW heat output of the modern Earth (and therefore the likely order of magnitude for non-tidal heat output from such world's silicate/iron cores even prior to tidal heating). On the other hand, eccentricity values well above e = 0.1 are observed in exosystems, and tidal heat (even for complex multilayer worlds) will scale as e². Similarly, many hosts smaller and less luminous than 55 Cnc exist, providing opportunities for ice mantles at shorter periods, and tidal heating strongly scales as the fifth power of mean motion, n⁵. Therefore, only a small increase in albedo, a small decrease in host luminosity, or simply the allowance of a thin liquid water surface ocean followed by clathrates and solid ices after only a few 100 km's depth can all lead to far greater heating. Lastly, Figure 8 uses a conservative ice viscosity of 1 × 10¹⁵, but lower values remain plausible due to issues discussed in Section 3.3.

The first two cases shown in Figure 8 illustrate an important point regarding Q values. For the same orbit, the first case with a 10,000 km ice mantle and Q_eff ∼ 9 produces nearly ten times the total heat in Watts compared to the second case with a 6000 km ice mantle and Q_eff ∼ 8. This occurs due to changes in the planetary Love numbers, and highlights the reason why the value −Im(k₂) for a given planet is a significantly better scalar measure of the world's tidal susceptibility than a Q value alone. Recall that −Im(k₂) replaces the ratio k₂/Q in the viscoelastic form of the classical tidal equation. For this reason, later figures here report −Im(k₂) values, which may be converted to effective Q values via: Q_eff = Re(k₂)/− Im(k₂). In later use, we generally omit the subscript eff, as all forms of Q values are already only a proxy for susceptibility to tidal heating, and only W(n, T) for a given planetary case is the final measure.

The Love and Shida numbers for simple ice–silicate worlds with ice mantle thicknesses up to 10,000 km are shown in Figure 10. Near and above this thickness, large gas envelopes are expected to alter results, unless unusual envelope stripping has occurred. Results are separated between the real and imaginary components of each Love number. In the Fourier model of viscoelastic tides, real components generally characterize elastic behavior, while imaginary components characterize energy dissipation. k₂, representing the gravitational response, is directly applied to tidal heat predictions. In a homogeneous case, −Im(k₂) may characterize the entire tidal heat response, and in these cases where the ice mantle dominates dissipation, it becomes a reasonable proxy for total tidal heating. The displacement Love number h₂ and Love–Shida number l₂ may be interpreted to signify the magnitude of radial and tangential surface displacements respectively. While l₂ is most properly referred to as the Shida number or Love–Shida number, for clarity below we do refer to the three related values k, h, and l by the common aggregate term Love numbers.

Zoom In Zoom Out Reset image size

Figure 10 demonstrates a material tidal resonance with ice layer depth, which typically occurs between t_ice = 1000–2000 km. Above these thicknesses, Love numbers grow or contract monotonically with increasing depth, with little change in magnitude for imaginary terms. Re(k₂) and Re(h₂) shift steadily upward for deep mantles but remain far from their hydrostatic limits of 3/2 and 5/2, respectively. Such limits are more typical of gravity dominated gas giants. Re(k₂) and Re(h₂) decrease relative to a 1 M_E core for ice depths below 2000 km, as the low density of ice allows the aggregate objects to briefly appear more strength dominated. Orbital period has a strong impact on the depth and especially the depth range that causes tidal resonance. Peak magnitudes in imaginary components remain similar across the whole period range shown. In particular, note how the full high resonant −Im(k₂) response may occur across the enormous ice thickness range from about 3000–10,000 km at short periods, while at long periods an ice mantle must be perfectly tuned by chance to around 1000 ± 100 km in order to experience a similarly high response.

While silicate-only mantles of this thickness are unrealistic, they are shown for comparison with ice results and represent a high-density, high-viscosity limit, which high pressure ice may approach. We plot constant density curves to help bound the behavior of the physical system. Full models of the interior from an equation of state and thermal flux model will, in addition to adding numerous new poorly constrained degrees of freedom to the system, be fully temperature dependent, and require iterative solutions with tidal heat flux as discussed in Section 3.3. As tidal heating depends far more strongly on viscosity than on density, the benefit of a full equation of state approach is also lessened when our goal is determining heat outcomes. Parameter gradations are discussed in more detail in Section 5.2, with results showing that low pressure upper ice levels are the primary control on dissipation.

While this paper does focus on questions of heating, very high surface displacement predictions are typical in strong tidal forcing cases, and displacement is significantly greater in low rigidity (4 × 10⁹ Pa) ice than in any underlying (5 × 10^{10 − 11} Pa) silicate. For example, lateral surface movements of ∼4 km, and radial displacements of ∼1 km may occur for a 4000 km thick ice mantle in a three-day period. This corresponds to surface strains of _θθ ∼ 0.005, _rr ∼ 0.003, _ϕϕ ∼ 0.002, and _θϕ ∼ 0.001. In general, the high strains experienced for such worlds suggest the possibly of extreme faulting, perhaps analogous to tidal features on Europa (Hurford et al. 2005, 2007b). The opening of faults by the tidal cycle may also induce periodic eruptions of volatiles as on Enceladus (Hurford et al. 2007a) and perhaps Europa (Roth et al. 2014), the timing of which may in the future be detectable through careful primary or secondary eclipse transit spectroscopy (Kaltenegger et al. 2010).

Fu et al. (2010) model the thermal structure of large ice mantles in detail. For the near surface of 2, 5, and 10 M_E super-Earths, they find a complex range of outcomes may arise from the interaction of the geotherm with the ice Ih, III, V, and VI portion of the phase diagram in the upper ∼60 km of planets. Both the presence or absence of a liquid ocean below a largely ice Ih crust is plausible, due to the unique notch in the liquid-solid phase boundary in the water phase diagram in this pressure regime at temperatures near 250–280 K, and is largely controlled by internal heat flux. Naturally, the surface temperature, itself a strong function of atmospheric state and solar flux, is the primary control on whether or not such an ocean extends up to the atmospheric interface as well. Since our current model assumes a solid surface, we focus on worlds where all water oceans exist below some amount of solid ice. These models are expected to also be close approximations of the total tidal dissipation and internal deformation of worlds with thin or modest liquid water oceans at the surface. In such an extended application, surface Love numbers here for the solid top case clearly do not represent the shape of a thin-to-modest liquid ocean's top surface that might be observed, but do approximate the shape of the base of such a layer.

For surface temperatures above 273 K, it is somewhat difficult to achieve a solid ice surface due to pressure alone. At 274 K, ∼6 kbar or 630 MPa is required, which would require an unrealistically large gas envelope for a 1.4 M_E world. Therefore, our analysis is limited to cases where a surface ocean, if present in cases of T_surf ⩾ 273 K, will be thin enough that Love numbers for the layer system beginning at the first solid interface are a close approximation of the body as a whole. 630 MPa is for example achieved under almost exactly 100 km of water for a 1.4 M_E world.

Beyond this, the criteria for a solid top layer is best achieved in cold environments. This is somewhat incompatible with strong tidal forcing of planets near to their host stars. Therefore, the analysis here is most useful for ice–silicate hybrid worlds encircling dim host stars, such as L and M-dwarfs, white dwarfs, neutron stars, or pulsars. Similarly, such a combination of strong tidal forcing and cold surface temperatures is achieved for larger moons of Jupiter and super-Jupiter class worlds, as well as the tidal evolution of binary terrestrial planet systems.

Water oceans may also exist at significant depths due to rising temperatures. The basal pressure of an ice mantle of ∼3000 km thickness, atop an Earth-mass core, is around 19 GPa, which for temperatures between ∼250 K to 750 K places ice in the ice VII phase. Above 750 K, supercritical liquid water is still possible at this depth (where all liquid above 647 K is in the supercritical state). The majority of a 3000 km ice mantle should participate in convection, except for worlds with unrealistically small silicate cores or worlds vastly older than the solar system timescales of interest here. The adiabatic geotherm in a convecting ice layer is typically shallow, making temperatures much beyond this range, even at depth, less likely. Eutectic melting point depression in ammonia-water systems has also been suggested as a common planetary mechanism for liquid ocean maintenance, and may allow for a complex diversity of liquid layers, including perched oceans with solid ices above and below.

The density of Ice VII in this deep and warm regime is still marginally well represented by ∼1000 kg m⁻³. Below ∼400 K the density climbs above 1100 kg m⁻³, but only reaches 1200 kg m⁻³ when as cold as around 200 K. Therefore, our focus may remain on viscosity variations and not density variations.

Figure 11 shows tidal surface Love numbers for ice–silicate hybrid super-Earth planets, showing variations as an ice mantle up to 3000 km thickness is added above a baseline Earth model of 1 R_E, 1 M_E, with a solid iron inner core, liquid outer core, and uniform silicate mantle. The maximum planet mass in these plots is 1.39 M_E, with a total radius of 1.47 R_E. A forcing period of P = 30 days is used. Note the log scales for imaginary components. In Case 11B a solid ice layer of density 1000 kg m⁻³, shear modulus 4 × 10⁹ Pa, and viscosity 1 × 10¹⁴ Pa s, is applied. This solid-only case leads to significant variations in h₂ and l₂ as a function of ice mantle thickness, and the highest dissipative response through −Im(k₂). In Case 11C this same ice mantle is separated from the silicate below by a 10 km thick water ocean. This greatly reduces the dissipative (imaginary) response, while h₂ and l₂ switch to monotonic behavior because they no longer reflect a transition from dominance by the silicate to dominance by the ice. Re(k₂), which depends largely on the density structure of layers, is minimally changed. In Case 11D this decoupling ocean is enlarged to 500 km. This has negligible impact on the dissipative response of the planet, and has only a minor impact on real valued components through the additional mass. Case 11E tests the concept of a perched water ocean, by placing a 10 km water layer at the midpoint of the given ice thickness. This significantly reduces the dissipative response, but has a real response nearly identical to Case 11C. Case 11F places a 10 km water ocean below a 20 km surface ice shell, with the bulk of the ice below both, directly attached to the silicate mantle. The result is similar to Case 11E but with less reduction in dissipation. For comparison, in Case 11A the silicate mantle is extended to match the overall radius of the ice-hybrid world, leading to a very low dissipative response, and a gradually more hydrostatic real-valued response. Overall this indicates that ocean presence is critical, ocean thickness negligible, and ocean position important mainly for the magnitude of dissipation.

Zoom In Zoom Out Reset image size

Figure 12 shows the effect of varying the orbital forcing period on the surface Love numbers of ice–silicate hybrid super-Earths. Baseline Earth model as in Figure 11. Solid curves: Ice mantles with no water ocean, varied from a 2 to 50 day orbital forcing period. At longer periods, the peak (Re(l₂), −Im(k₂)), or minimum (Re(k₂), Re(h₂)) response of non-monotonic functions occurs at smaller ice mantle thicknesses. Dashed Curves: Ice mantles are separated from the silicate mantle below by a 10 km water ocean. Changes in real-valued components are negligible with forcing frequency, while the dissipative response (−Im(k₂)) is highest at short periods. This is the opposite of the dissipative trend for the solid cases for thinner ice mantles (prior to any peak), which achieve maximum dissipation at longer periods. Note that for decoupled ice mantles, the shape of the imaginary component curves is entirely insensitive to the forcing period, and thus dissipation is always greater with greater thickness, while for solid-only ice mantles, there is generally a thickness for peak dissipation at any given forcing frequency.

Zoom In Zoom Out Reset image size

Continuous variation of ocean position was tested and found to have minor impact beyond that already reflected by base and near surface ocean position endmembers. There is some dependence of Re(k₂) and Re(h₂) on position, but it is weak and mainly expressed at short periods, such as P = 2 days (Note that Love numbers are a function of the planet only, and therefore apply to all host masses and eccentricities). At P = 2 days, all ocean positions lead to increasing −Im(k₂) with increasing ice thickness. At longer periods, near surface oceans lead to decreasing −Im(k₂) with ice thickness. At long periods, e.g., P = 300 days, even the dependence of −Im(k₂) on ocean position begins to be lost, but the dependence of Re(l₂) is enhanced. Near surface oceans lead to the highest ocean bearing branch for Re(l₂), consistent with the idea that freedom of an outer ice shell to deform laterally without connection to the lower solid increases the overall lateral crustal deformation response.

When a uniform growing ice mantle is applied to a silicate-iron core with gradations in the silicate properties using the PREM baseline, results almost exactly overlie those for a single or double layer silicate mantle model as shown in Figure 11, with the exception of slight constant offsets due to any imperfectly matched ice-free values. We find this to be true both at short and long forcing periods. Thus, the importance of such iron/silicate interior gradations is negligible in terms of surface Love numbers for ice hybrid worlds, except when ice layers are less than a few hundred km in thickness.

Of greater importance is any gradation of ice viscosities. Cold brittle ice with a viscosity up to 1 × 10²¹ Pa s has been suggested for the non-convective top layer of Europa's ice shell, yet near surface warm convective ice may have viscosities from 1 × 10¹³ to 1 × 10¹⁵ Pa s depending on composition and conditions. As with silicates, there are competing effects upon the viscosity of ice with depth: as high pressure may act to increase or decrease viscosities, while increasing temperature with depth acts to lower them. Thus, it is not immediately clear for a given composition and global heat flux which behavior may dominate, as discussed in Section 3.3.

Compared to the range of variation possible for ice viscosities, the range available for ice densities and ice shear moduli, even under high pressures in a deep ice mantle, are comparatively small. Ice density may increase by a few percent, but it will not vary over up to two orders of magnitude. Therefore Figure 13 shows models where only the viscosity is varied in the ice layer, so as to deconvolve this strong primary effect from any other more modest parameter variations.

Zoom In Zoom Out Reset image size

In general, ice property gradations have a major impact on outcomes, however, it is usually the uppermost layer properties which continue to dominate the response. Because the overall viscosity has changed, imaginary components experience significant changes in magnitude. A brittle top suppresses dissipation (e.g., lowers −Im(k₂) values) both with and without an ocean (13D and 13H). Unlike for solid-only cases, for models with a thin ocean, both forms of a viscosity gradation (13F increasing with depth from 1 × 10¹³ to 1 × 10¹⁵, or 13G decreasing with depth from 1 × 10¹⁵ to 1 × 10¹³), will increase dissipation. Here amplification due to the lowermost layer is also activated partly due to the higher stress concentration typical at interfaces between a solid top layer and underlying ocean, as seen in Figure 4. Gradations also alter system natural frequencies, such that the ice thickness for peak dissipation or maximal or minimal deformation shifts from a pure case.

The trend for ice gradations to have high importance continues at other forcing periods. Results for a short period of three days are shown in Figure 14, and are distinctly different from patterns in Figure 13. At this short period a much greater fraction of structure cases lead to dissipation which increases monotonically with ice thickness, and overall lowers the chances that any response peak of a property lies within the 0–3000 km ice mantle range.

Zoom In Zoom Out Reset image size

The most interesting gradation case occurs at long periods, where a double peak in the dissipative response appears in the case of a solid ice mantle with a brittle top. In this scenario the natural frequency of the convective ice layer and brittle ice layer are sufficiently unique to produce two visibly separate peaks, and the forcing period of 100 days in Figure 13 happens to excite both within the ice thickness range applied. This example is illustrative of a general phenomenon that may occur for high viscosity contrast layers. A similar tendency to multiple peaks does not occur for both the silicate and ice layers due to the challenge of exciting a thick deep 1 × 10²¹ Pa s layer versus the easier task of exciting such a brittle layer as a thin surface shell.

Figure 15 concludes by showing the tidal behavior for selected models as functions of orbital period. Total heat, effective Q values, and circularization times following Equation (2) are shown. Again a host star mass and luminosity following 55 Cnc is used to allow colder surfaces so that direct comparisons may be made between ice and silicate worlds, however, curves in Figure 15 are nearly identical to those for a host star with a mass of 1 M_sol. Silicate-iron worlds are shown in the top row, while ice–silicate worlds are shown in the bottom row. Cases are selected to highlight the role of low viscosities, magma oceans, and water oceans. Even with a conservative assumption on eccentricity, geologically relevant tidal heating is shown to extend out to ∼50–100 days in a variety of cases, with a horizontal threshold of 5 TW used to denote this approximate activity level. Geologically relevant tidal activity extends to longer periods more often in ice–silicate hybrid cases (bottom left panel), with deep-ice mantle cases (e.g., 6000 km of ice or more) constituting the longest period active worlds.

Zoom In Zoom Out Reset image size

For Earth-mass silicate-iron worlds, all structure variations tested in Figure 4 cluster close to the second curve in Figure 15 labeled 1 M_E Multilayer, including both PREM variants. The most important result of Figure 15 is the change that occurs between worlds we label as Cool Dry Earths, Warm Earths, and Very Hot Earths. By Cool Dry Earths we refer to worlds without dissipation in a surface water ocean, with mantle structures and mantle viscosities similar to the modern Earth. The high viscosity of such mantles leads to very low dissipation, very large Q values, and very long circularization times. Such worlds may retain initial eccentricities for much longer than solar system lifetimes, even at short periods, given that their eccentricities are not so large as to heat them into the next Warm Earth category. This finding may help to explain some high eccentricities observed in the exoplanet population for terrestrial class planets.

A significant change occurs for Warm Earth cases, when we assume the mantle and asthenosphere viscosities of silicate material are far lower than for the modern Earth, and close to viscosities at or just above typical solidus temperatures for mantle materials. For a mantle viscosity of 1 × 10¹⁷ Pa s (e.g., high pressure, just sub-solidus) and asthenosphere viscosity of 1 × 10¹⁴ Pa s (e.g., just above solidus), tidal heating of ∼5 TW may occur out to 80 day periods at e = 0.1, to 320 days for e = 0.2, and 2000 days for e = 0.5. Solutions of ∼10,000 TW are plausible for partially melted silicate worlds out to periods of 20, 80 and 500 days for the same eccentricities. While higher order terms in eccentricity become important above this level, such moderate Warm Earth silicate worlds have the potential to maintain the heat required for a high partial melt fraction mantle across a very large range of orbital periods relevant for terrestrial planet formation. Following Figure 1, the definition of Warm Earth above corresponds to bulk layer temperatures that are only ∼300–400 K higher than the present Earth. Effective Q values plummet for Warm Earths, falling 4–5 orders of magnitude relative to a Cool Dry Earth, and well below the common assumption of Q = 100. Typical Warm Earth Q values are in the range 1–10, with values of Q = 1–2 very common for short periods. Circularization times for such worlds are far faster than their cool counterparts, a result that can help warm Earth-analog worlds to rapidly recover from high eccentricity scattering events, and reduce the number of terrestrial planets vulnerable to being lost from early active solar systems due to orbit crossings or three body encounters.

Only for the most extreme Very Hot Earths, does the opposite behavior occur, with Q values again rising back above 100, and into the Q = 1000 range. We emphasize that this is a highly unexpected result. Conventional wisdom suggests that melting in general reduces dissipation, thus extending circularization times, and making Earths vulnerable to orbital crossings and system loss for longer times. Instead, we find that the path of multilayer global heating first produces an extreme increase in dissipation, thus shorter circularization times, and only once a magma ocean consumes a very large fraction of a mantle, does traditional tidal shutdown begin. This is a uniquely multilayer effect, caused by the low viscosity warm mantle underneath a magma ocean. Here, even the imperfect viscoelastic tuning of warm lower mantle material more than compensates for the loss of active volume from magma layers. Note how even the most extreme magma ocean cases shown have lower Q values (greater dissipation) than Cool Dry Earths, again counter to the notion of tidal shutdown. Overall, all dry planet results with Q near the 20–100 level are rare, and thus the assumption of Q ∼ 20–100 is most useful for Cool Wet Earths: e.g., planets expected to strongly resemble the modern Earth in age, orbital forcing, and expected surface water ocean state.

The viscosity at which behavioral transitions occur is period dependent. A Maxwell time of 100 days is achieved for silicate material with a viscosity of around 1 × 10¹⁷ Pa s (or 1 × 10¹⁶ Pa s for 10 days). At this period, below this viscosity, the tidal response is again attenuated. For the cases containing magma oceans in Figure 15 we have tested values for the remnant solid lower mantle of 1 × 10¹²–1 × 10¹⁴ Pa s, to try to bring about the least dissipation plausible for such layers while remaining a solid. Lower mantle viscosities any closer to 1 × 10¹⁶–1 × 10¹⁷ Pa s in Very Hot Earth cases will show even less of an increase in Q. The fact that creating Very Hot Earth class results requires stretching the parameter range to this extent hints at how common Warm Earth class results will be. The relative position of Q values for each category can be viewed as an aggregate reflection of the distance from the ideal Maxwell tuning of each worlds' dominant layers: e.g., at ∼1 × 10¹⁷ Pa s, Warm Earths are ideally tuned, while Very Hot Earths at ∼1 × 10¹⁴ Pa s are still typically closer to their ideal tuning than Cool Dry Earths at ∼1 × 10²¹ Pa s. For improved rheologies beyond the Maxwell model, specific ideal tuning frequencies will vary, often with more gradual transitions between categories, but regardless of rheology the overall pattern will remain robust that dissipation is lowest for cool and hot cases, and highest for intermediate cases.

Previous analysis (Henning et al. 2009) has shown that feedback with mantle convection by the mechanism described in Moore (2003b) will cause the majority of strong tidal forcing cases to rapidly evolve into very stable tidal-convective equilibrium states at low-to-modest melt fractions. While this previous work was performed for a homogeneous mantle, the general tidal-convective equilibrium mechanism makes it very difficult for planets to reach the Very Hot Earth states with deep magma oceans demonstrated here. We conclude such states are only possible when external orbital forcing is extremely strong, such as for P ⩽ 10 days, or at longer periods e ⩾ 0.5. This means that the vast majority of terrestrial planets warmed any significant amount more than the modern Earth, are expected to fall within the Warm Earth category, and therefore will have circularization times 10–100 times faster than Q = 100 models would predict.

Transitions between our categories of Cool Dry Earth, Warm Earth, and Very Hot Earth are gradual and difficult to define, based on the large number of parameters involved, and sensitivity to forcing frequency in multiple layers. We have highlighted approximate bounding cases, to collapse parametric variants into a single behavioral spectrum, roughly analogous to bulk planetary temperature. The temperature of a given exoplanet will be influenced by its age, orbital state, orbital history, internal structure, composition, and internal structure history. Eccentricity and semi-major axis however do provide useful guides, as low eccentricity and longer period objects are far more likely to remain cool, while the forcing to create Very Hot cases will arise primarily in high eccentricity and short-period objects. Note that the control due to semi-major axis arises from the n⁵ component of Equation (1), and not necessarily from a rise in isolation driven surface temperatures for short period worlds, as mantle temperatures are more controlled by total heat flux through convection than by the surface temperature boundary condition.

For ice–silicate hybrid worlds, Q ⩽ 100 for the majority of models, and is often ⩽10. Because of its low viscosity, the Maxwell time for ice is often crossed in the period range shown, leading to Q(P) curves with distinct minima that shift based on model details. Overall however, there is no major difference in behaviors when small water oceans are added to previously all-solid models. The main way to achieve high dissipation of ice systems is to increase ice mantle depth. The primary pathway to decrease dissipation is to remove ice volume by converting it to liquid. Other variants, such in viscosity or ocean location, continue to cluster near the baseline model shown here with no simple relationship. The importance of ice for tidal modeling is therefore the basic extension of tidal relevance to long periods (100 days) for colder worlds, compared to the very short period range of tidal relevance for cold silicate worlds (10 days). The inclusion of water oceans does not alter this result.

We find that modest global-scale partial melting in effect causes silicate and ice worlds to switch tidal behaviors. Normally low-dissipative silicate worlds become high-dissipative, and normally high-dissipative ice worlds become low-dissipative. Full role reversal, however, is more common at long periods, and extreme melting (e.g., magma or water oceans spanning ∼1/2 or more of a mantle layer's volume, perhaps driven by extreme orbits or concurrent spin-tides) in all worlds causes low-dissipative behavior.

The low-Q high-dissipation behavior of super-Earths with large ice mantles implies that a major transformation must occur when worlds grow to the size of Neptune where Q = 12,000–330,000 (Banfield & Murray 1992). This transition may be dominated by very high pressure ice phases such as superionic ice and plasma ice (Ryzhkin 1985; Goldman et al. 2005; French et al. 2009; Goncharov et al. 2009) expected to exist in Neptune's mantle (D. Sasselov 2011, private communication). Recent research (Cavazzoni et al. 1999; Redmer et al. 2011) suggests that Neptune's geotherm may pass very close to the transition between the cooler superionic phase and the warmer plasma phase, possibly dividing the planet's mantle into both materials. The viscosity of these phases is unknown, but the characteristic that superionic ice has oxygen nuclei retained in a lattice, while hydrogen nuclei are mobile, yet in plasma ice both species become mobilized, suggests that the warmer plasma phase may have viscosities more like a fluid. Therefore as Neptune-class worlds cool in time, their mantles may switch from a warmer fluid-like tidal behavior with high Q values, to a possible evolved era at lower Q when more superionic ice forms, potentially akin to the 10,000 km thick ice mantle models in this work, yet depending on future determinations of high pressure viscosities.