Tidal Heating of Earth-like Exoplanets around M Stars: Thermal,Magnetic,and Orbital Evolutions

Abstract

The internal thermal and magnetic evolution of rocky exoplanets is critical to their habitability. We focus on the thermal-orbital evolution of Earth-mass planets around low-mass M stars whose radiative habitable zone overlaps with the “tidal zone,” where tidal dissipation is expected to be a significant heat source in the interior. We develop a thermal-orbital evolution model calibrated to Earth that couples tidal dissipation, with a temperature-dependent Maxwell rheology, to orbital circularization and migration. We illustrate thermal-orbital steady states where surface heat flow is balanced by tidal dissipation and cooling can be stalled for billions of years until circularization occurs. Orbital energy dissipated as tidal heat in the interior drives both inward migration and circularization, with a circularization time that is inversely proportional to the dissipation rate. We identify a peak in the internal dissipation rate as the mantle passes through a viscoelastic state at mantle temperatures near 1800 K. Planets orbiting a 0.1 solar-mass star within 0.07 AU circularize before 10 Gyr, independent of initial eccentricity. Once circular, these planets cool monotonically and maintain dynamos similar to that of Earth. Planets forced into eccentric orbits can experience a super-cooling of the core and rapid core solidification, inhibiting dynamo action for planets in the habitable zone. We find that tidal heating is insignificant in the habitable zone around 0.45 (or larger) solar-mass stars because tidal dissipation is a stronger function of orbital distance than stellar mass, and the habitable zone is farther from larger stars. Suppression of the planetary magnetic field exposes the atmosphere to stellar wind erosion and the surface to harmful radiation. In addition to weak magnetic fields, massive melt eruption rates and prolonged magma oceans may render eccentric planets in the habitable zone of low-mass stars inhospitable for life. Key Words: Tidal dissipation—Thermal history—Planetary interiors—Magnetic field. Astrobiology 15, 739–760.

1. Introduction

Gravitational tides are common in the Solar System, from the Moon, responsible for driving the principle diurnal tides in Earth's oceans and atmosphere, to Io, the most volcanically active body in the Solar System. Tidal dissipation as a heat source in the solid Earth is weak at present and often neglected from thermal history calculations of its interior. However, rocky exoplanets with eccentric orbits close to their star are expected to experience significant tides (Dole, 1964; Rasio et al., 1996; Jackson et al., 2009; Barnes et al., 2010) that likely influence their thermal, orbital, and even atmospheric evolution (Barnes et al., 2013; Luger et al., 2015). Recent advances in modeling tidal dissipation in a viscoelastic mantle (Henning et al., 2009; Běhounková et al., 2010, 2011; Henning and Hurford, 2014) have advocated using Maxwell-type temperature-pressure-dependent rheology and emphasized the limited applicability of a constant tidal quality factor “ \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${\mathscr Q}$$ \end{document} ” model. These more complicated dissipation models are necessary to better characterize the tidal and orbital states of rocky exoplanets over a range of internal temperatures.

The search for habitable Earth-like exoplanets commonly targets planets in orbit around low-mass M type stars to maximize the number of small-mass planets found (Seager, 2013; Mayor et al., 2014). Targeting low-mass stars is beneficial for at least three reasons: (1) the habitable zone (HZ) around M stars is much closer to the star (Kopparapu et al., 2013), making an Earth-mass planet in the HZ an easier target for both transit and radial velocity detection, (2) low-mass M stars are more abundant in the nearby solar neighborhood, and (3) M stars have longer main sequence times.

On the other hand, there are several reasons why targeting M stars may be risky. M stars are intrinsically faint, which makes most observations low signal to noise, and their flux peaks at wavelengths close to those absorbed by Earth's atmosphere. M stars are more active, especially early on, which may induce massive amounts of atmospheric loss (Luger et al., 2015) and biologically hazardous levels of radiation at the surface. Earth-mass planets in the HZs of M stars likely experience larger gravitational tides associated with star-planet and planet-planet interactions, especially considering that most exoplanet systems are dynamically full (Barnes and Quinn, 2004; Barnes and Greenberg, 2006; Van Laerhoven et al., 2014). However, the implications of these tides on the thermal evolution of the interior have not yet been explored.

The thermal history of a planet is critical to its habitability. Mantle temperature determines the rates of melting, degassing, and tectonics, while the thermal state of the core is critical to the maintenance of a planetary magnetic field that shields the surface from high-energy radiation. The thermal evolution of Earth and terrestrial planets involves solving the time evolution of mantle and core temperature through a balance of heat sources and sinks. The thermal history of Earth, although better constrained than any other planet, is still subject to significant uncertainties. However, avoiding both the thermal catastrophe in the mantle (Korenaga, 2006) and the new core paradox (Olson et al., 2013) adds significant constraints that predict a monotonic cooling of the mantle and an active geodynamo over the history of the planet (Driscoll and Bercovici, 2014). Previous models of the thermal evolution of rocky exoplanets (e.g., Gaidos et al., 2010; Driscoll and Olson, 2011; Tachinami et al., 2011; Van Summeren et al., 2013; Zuluaga et al., 2013) have focused on the influence of planet size on thermal evolution but neglected tidal dissipation as an internal heat source; therefore magnetic field strength and lifetime in the HZ around low-mass stars were likely overestimated. Here we improve the thermal and magnetic evolution model of Driscoll and Bercovici (2014) by adding tidal heating as an internal heat source, and couple this to orbital evolution.

In this paper, we focus on the influence of tidal dissipation in the solid mantle of Earth-like exoplanets in the HZ around M stars. Tidal dissipation deposits heat in the planetary interior, while simultaneously extracting energy from the orbit, which can lead to circularization and migration. We couple the thermal and orbital evolution equations into a single model to identify the conditions and timescales for Earth-like geophysical and magnetic evolution. Section 2 describes the thermal-orbital model equations. Steady-state behavior is discussed in Section 3 to build intuition about the thermal-orbital coupling. Results with evolving orbits are presented in Section 4 and with fixed orbits, mimicking forcing by companion planets, in Section 5. The possibility of an internally driven runaway greenhouse is addressed in Section 6. The influence of tides on the inner edge of the HZ over a range of stellar masses is explored in Section 7. A summary and discussion are in Section 8.

2. Model Description

In this section, we describe the details of the thermal-orbital evolution model. Section 2.1 describes the tidal dissipation model, Section 2.2 describes the thermal evolution of the coupled mantle-core interior, and Section 2.3 describes the orbital evolution as a function of dissipation efficiency. A note on terminology: here Q's refer to heat flows (units of W), q's refer to heat fluxes (units of W m⁻²), and the script symbol \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\cal{Q}$$ \end{document} refers to the tidal quality factor (dimensionless), also known as the specific dissipation function.

2.1. Tidal dissipation model

The gravitational perturbation experienced by a secondary body (the planet) in orbit about a primary body (the star) is approximated by the lowest-order term in the potential expansion, which is the semidiurnal tide of degree 2 (Kaula, 1964). The power dissipated by tidal strain associated with this term in the secondary with synchronous rotation is (Segatz et al., 1988) \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}Q_ { \rm tidal } = - \frac { 21 } { 2 } { \rm Im } ( k_2 ) \frac { GM_*^2 { R_ { \rm p } } ^5 \omega e^2 } { a^6 } \tag { 1 } \end{align*} \end{document}

where G is the gravitational constant, M _* is stellar mass, R _p is planet radius, ω is orbital frequency, e is orbital eccentricity, a is orbital semimajor axis, and Im(k ₂) is the imaginary part of the complex second-order Love number k ₂. If planetary rotation is synchronous, then the tidal frequency is equal to the mean motion \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\omega = n = \sqrt{GM_* / a^3}$$ \end{document} , and the tidal power reduces to \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}Q_ { \rm tidal } = - \frac { 21 } { 2 } { \rm Im } ( k_2 ) G^ { 3 / 2 } M_*^ { 5 / 2 } { R_ { \rm p } } ^5 \frac { e^2 } { a^ { 15 / 2 } } \tag { 2 } \end{align*} \end{document}

This expression for tidal dissipation is the product of three physical components: (1) tidal efficiency (−Im(k ₂)), (2) star-planet size \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$( M_*^{5 / 2} {R_{\rm p}}^5 )$$ \end{document} , and (3) orbit (e ²/a ^15/2).

For illustrative purposes, it is helpful to compare the radiative “habitable” zone (HZ) to the “tidal” zone, the orbital distance at which tidal dissipation is likely to dominate the internal heat budget of the planet. For this comparison, we compute the tidal heat flow using (2) for an Earth-sized planet with e = 0.1 and −Im(k ₂) = 3 × 10⁻³, similar to present-day Earth. Figure 1 shows that for M _star < 0.3 M _* the HZ overlaps with the tidal zone, defined as when tidal heating is the dominant term in the heat budget. This diagram implies that Earth-mass planets in the HZ of low-mass stars could experience extreme tidal heating and a rapid resurfacing rate that may render the surface uninhabitable. The HZ and tidal zone intersect because stellar radiation flux is more sensitive to stellar mass than the gravitational tidal potential. We note that larger eccentricity or tidal dissipation efficiency (−Im(k ₂)) would push the tidal zone limits out to larger orbital distances, rendering the HZ tidally dominated for larger-mass stars.

FIG. 1.

Comparison of the radiative “habitable” zone to the “tidal” zone. The radiative “habitable” zone is from Kopparapu et al. (2013). Inside the “tidal” zone, heat released by tidal dissipation is likely to dominate the internal heat budget of the planet. The tidal zone is delineated by distances from the star where an Earth-mass planet would receive an amount of heat via tidal dissipation equal to either the surface heat flow of Io (Q _surf ≈ 80 TW, left curve) or Earth (Q _surf = 40 TW, right curve). Tidal heat flow is calculated by (2) assuming e = 0.1 and −Im(k ₂) = 3 × 10⁻³ (k ₂ = 0.3, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\cal{Q}$$ \end{document} = 100). The gray shaded region denotes the zone where the planet is predicted to be radiatively “habitable” but tidally dominated, and therefore possibly not habitable.

The one-dimensional dissipation model in (2) assumes a homogeneous body with uniform stiffness and viscosity. To derive the dissipation efficiency (−Im(k ₂)), we first define the Love number, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}k_2 = \frac {3 / 2 } {1 + \frac {19} {2} \frac {\mu}{\beta_{st}}} \tag{3}\end{align*} \end{document}

where μ is shear modulus and β _st is effective stiffness (Peale and Cassen, 1978). Writing shear modulus as a complex number and using the constitutive relation for a Maxwell body (Henning et al., 2009), one can derive the dissipation efficiency in (2) as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}- { \rm Im } ( k_2 ) = \frac { 57 \eta \omega } { 4 \beta_ { \rm st } \left( 1 + \left[ \left( 1 + \frac { 19 \mu } { 2 \beta_ { \rm st } } \right) \frac { \eta \omega } { \mu } \right] ^2 \right) } \tag { 4 } \end{align*} \end{document}

where η is dynamic viscosity. We note that this model does not involve a tidal \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\cal {Q}$$ \end{document} factor; rather, the rheological response of the mantle is described entirely by Im(k ₂). For comparison with other models, one can compute the standard tidal \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\cal {Q}$$ \end{document} factor of the Maxwell model as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\cal {Q} = \frac { \eta \omega } { \mu } \tag { 5 } \end{align*} \end{document}

The common approximation is then −Im(k ₂) ≈ Re(k ₂)/ \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\cal {Q}$$ \end{document} (e.g., Goldreich and Soter, 1966; Jackson et al., 2009; Barnes et al., 2013). Although the tidal \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\cal {Q}$$ \end{document} factor does not appear explicitly in the calculations below, it will be used to calibrate this model to Earth and is useful in comparing it to other models (also see Appendix A).

Following previous studies (Sotin et al., 2009; Běhounková et al., 2010, 2011), we ensure that the model reproduces the observed tidal dissipation in the solid Earth by calibrating the effective material properties in (4) appropriately. This calibration allows us to approximate the total tidal dissipation over the whole mantle by a single volume-averaged dissipation function. The \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\cal {Q}$$ \end{document} factor of the solid Earth has been estimated empirically to be \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\cal {Q}$$ \end{document} _E≈100 (Ray and Egbert, 2012). Effective viscosity follows an Arrhenius law form, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\nu = \nu_ { \rm ref } { \rm exp } \left( \frac { E_ { \nu } } { R_ { \rm g } T_ { \rm m } } \right) / \epsilon_ { \rm phase } \tag { 6 } \end{align*} \end{document}

where ν = η/ρ _m is kinematic viscosity, ρ _m is mantle density, ν _ref is a reference viscosity, E_ν is the viscosity activation energy, R _g is the gas constant, T _m is average mantle temperature, and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\epsilon_{ \rm phase}$$ \end{document} accounts for the effect of the solid to liquid phase change (see Table 1 for a list of constants). Shear modulus is similarly described, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\mu = \mu_ { \rm ref } { \rm exp } \left( \frac { E_ { \mu } } { R_ { \rm g } T_ { \rm m } } \right) / \epsilon_ { \rm phase } \tag { 7 } \end{align*} \end{document}

Table 1.

Model Constants

Symbol	Value	Units	Reference
A_ν	3 × 10⁵	J mol⁻¹	Viscosity activation energy in (6)
A_μ	2 × 10⁵	J mol⁻¹	Nominal shear modulus activation energy in (7)
A _sol	−1.160 × 10⁻¹⁶	K m⁻³	Solidus coefficient in (29) (ET08)
α	3 × 10⁻⁵	K⁻¹	Thermal expansivity of mantle
α _c	1 × 10⁻⁵	K⁻¹	Thermal expansivity of core
B	2.5	nd	Melt fraction coefficient in (8)
B _sol	1.708 × 10⁻⁹	K m⁻²	Solidus coefficient in (29), calibrated
β	1/3	nd	Convective cooling exponent in (25)
β _st	1.71 × 10⁴	GPa	Effective mantle stiffness, calibrated in Section 2.1
c _m	1,265	J kg⁻¹ K⁻¹	Specific heat of mantle
c _c	840	J kg⁻¹ K⁻¹	Specific heat of core
C _sol	−9.074 × 10⁻³	K m⁻¹	Solidus coefficient in (29), calibrated
D	2,891	km	Mantle depth
D _Fe	7,000	km	Iron solidus length scale
D_N	6,340	km	Core adiabatic length scale
D _sol	1.993 × 10⁴	K	Solidus coefficient in (29), calibrated
δ _ph	6	nd	Rheology phase coefficient in (8, 9)
E _G	3 × 10⁵	J kg⁻¹	Gravitational energy density release at the ICB
\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\epsilon$$ \end{document} _UM	0.7	nd	Upper mantle adiabatic temperature drop
\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\epsilon$$ \end{document} _LM	1.3	nd	Lower mantle adiabatic temperature jump
\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\epsilon$$ \end{document} _c	0.8	nd	Average core to CMB adiabatic temperature drop
φ ^*	0.8	nd	Rheology phase coefficient in (8, 9)
g _UM	9.8	m s⁻²	Upper mantle gravity
g _LM	10.5	m s⁻²	Lower mantle gravity
g _c	10.5	m s⁻²	CMB gravity
γ _core	1.3	nd	Core Gruneisen parameter
γ _dip	0.2	nd	Magnetic dipole intensity coefficient in (35)
γ _ph	6	nd	Rheology phase coefficient in (8, 9)
k _UM	4.2	W m⁻¹ K⁻¹	Upper mantle thermal conductivity
k _LM	10	W m⁻¹ K⁻¹	Lower mantle thermal conductivity
κ	10⁶	m² s⁻¹	Mantle thermal diffusivity
L _Fe	750	kJ kg⁻¹	Latent heat of inner core crystallization
L _melt	320	kJ kg⁻¹	Latent heat of mantle melting
L _c	2.5 × 10⁻⁸	WΩK⁻¹	Lorentz number
L _^*	3.09 × 10²³	W	Stellar luminosity for M _^* = 0.1M _sun (B13)
M _m	4.06 × 10²⁴	kg	Mantle mass
M _c	1.95 × 10²⁴	kg	Core mass
μ _ref	6.24 × 10⁴	Pa	Reference shear modulus in (7)
μ ₀	4π × 10⁻⁷	H m⁻¹	Magnetic permeability
ν ₀	6 × 10⁷	m² s⁻¹	Reference viscosity
ν _LM /ν _UM	2	nd	Viscosity jump from upper to lower mantle
Q _rad,0	60	TW	Initial mantle radiogenic heat flow (J07)
R	6,371	km	Surface radius
R _c	3,480	km	Core radius
R _m	4,925	km	Radius to average mantle temperature T _m
Ra _c	660	nd	Critical Rayleigh number
ρ _c	11,900	kg m⁻³	Core density
ρ _ic	13,000	kg m⁻³	Inner core density
ρ _m	4,800	kg m⁻³	Mantle density
ρ _melt	2,700	kg m⁻³	Mantle melt density
ρ _solid	3,300	kg m⁻³	Mantle upwelling solid density
Δ ρ_χ	700	kg m⁻³	Outer core compositional density difference
σ _c	10 × 10⁵	S m⁻¹	Core electrical conductivity
T _Fe,0	5,600	K	Iron solidus coefficient in (23)
τ _rad	2.94	Gyr	Mantle radioactive decay timescale
τ _rad,c	1.2	Gyr	Core radioactive decay timescale
ξ	5 × 10⁻⁴	nd	Rheology phase coefficient in (8, 9)

Nondimensional units are denoted nd. References are B13 = Barnes et al. (2013); ET08 = Elkins-Tanton (2008); J07 = Jaupart et al. (2007).

This model predicts the rapid drop in shear modulus with melt fraction demonstrated experimentally by Jackson et al. (2004). The reference shear modulus μ _ref = 6.24 × 10⁴ Pa and effective stiffness β _st = 1.71 × 10⁴ GPa are calibrated by k ₂ = 0.3 and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\cal {Q}$$ \end{document} = 100 for the present-day mantle.

We model the influence of melt fraction φ on viscosity following the parameterization of Costa et al. (2009), \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\epsilon_ { \rm phase } ( \phi ) = \frac { 1 + \Phi^ { { \delta } _ { ph } } } { [ 1 + F ] ^ { B \phi^* } } \tag { 8 } \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}F = ( 1 - \xi ) { \rm erf } \left[ \frac { \sqrt { \pi } } { 2 ( 1 - \xi ) } \Phi ( 1 + \Phi^ { { \gamma } _ { ph }} ) \right] \tag { 9 } \end{align*} \end{document}

where Φ = φ/φ* and φ*, ξ, γ_ph , and δ_ph are empirical constants (Table 1).

The functional form of tidal shear modulus at high temperature-pressure is not well known, so we investigate the influence of shear modulus activation energy A_μ on the model. Contoured in Fig. 2 are tidal power using (2–4), tidal power using the approximation −Im(k ₂) = k ₂/ \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\cal {Q}$$ \end{document} , Love number k ₂, and tidal factor \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\cal {Q}$$ \end{document} as functions of shear modulus μ and viscosity ν. Evolution paths of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\cal {Q}_{ \rm tidal}$$ \end{document} as a function of T _m in the range 1500–2000 K are also shown in these contours for three shear modulus activation energies: A_μ = 0, 2 × 10⁵, and 4 × 10⁵ J mol⁻¹. For the nominal activation energy of A_μ = 2 × 10⁵ J mol⁻¹, the mantle cooling path passes through a maximum tidal dissipation around ν = 10¹⁶ Pa s and μ = 10¹⁰ Pa s, corresponding to T _m ≈ 1800 K where the mantle is in a viscoelastic state. In this state, the short-term tidal response of the mantle is elastic, while the long-term response is viscous (see Section 3). We note that, although these paths pass through a local maximum in \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\cal {Q}_{ \rm tidal}$$ \end{document} , they do not pass through the global maximum (dark red), which has been invoked for Io (Segatz et al., 1988). The main influence of increasing A_μ is to shift the dissipation peak to lower temperatures. The nominal value of A_μ = 2 × 10⁵ J mol⁻¹ produces a dissipation peak when melt fraction is about 50%.

FIG. 2.

Comparison of tidal properties as a function of viscosity and shear modulus for three shear-modulus activation energies A_μ = 0 J mol⁻¹ (blue line), 2 × 10⁵ J mol⁻¹ (black line), and 4 × 10⁵ J mol⁻¹ (red line). Lines show tracks of ν(T _m) and μ(T _m) for mantle temperatures in the range 1500–2000 K. (a) Contour of tidal heat flow Q _tidal. (b) Contour of tidal power using the approximation −Im(k ₂) = k ₂/ \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\cal {Q}$$ \end{document} . (c) Contour of Love number k ₂. (d) Contour of tidal dissipation factor \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\cal {Q}$$ \end{document} . Calculations use M _* = 0.1M _sun, A_ν = 3 × 10⁵ J mol⁻¹, e = 0.1.

The tidal dissipation factor \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\cal {Q}$$ \end{document} (Fig. 2d) and tidal power when using \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\cal {Q}$$ \end{document} (Fig. 2b) change by about one order of magnitude over the entire temperature and activation energy range (Fig. 2d). However, tidal power when using the Maxwell model in (2) fluctuates by 10⁵ TW over the same range (Fig. 2a). This comparison emphasizes that, although \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\cal {Q}$$ \end{document} is not far from constant, dissipation when using the Maxwell model is significantly different than the \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\cal {Q}$$ \end{document} approximation (also see Appendix A). The Love number k ₂ increases monotonically with T _m for all cases (Fig. 2c) up to the limit of k ₂ = 3/2 when μ/β << 2/19 in (3).

We note that dissipation in this model is a lower bound as dissipation in the liquid is not included, which can occur by resonant dissipation (e.g., Matsuyama, 2014; Tyler, 2014). Dissipation in the liquid is not likely to be a major heat source but could drive mechanical flows in the core (Zimmerman et al., 2014; Le Bars et al., 2015) and amplify dynamo action there (Dwyer et al., 2011; McWilliams, 2012).

2.2. Thermal evolution model

The thermal evolution of the interior solves the balance of heat sources and sinks in the mantle and core. The thermal evolution is modeled as that of Driscoll and Bercovici (2014), with an updated mantle solidus and inclusion of latent heat release due to magma ocean solidification (see Appendix B3). The conservation of energy in the mantle is \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}Q_{\rm surf} = Q_{\rm conv} + Q_{\rm melt} = Q_{\rm rad} + Q_{\rm cmb} + Q_{\rm man} + Q_{\rm tidal} + Q_{\rm L , man} \tag{10}\end{align*} \end{document}

where Q _surf is the total mantle surface heat flow (in W), Q _conv is heat conducted through the lithospheric thermal boundary layer that is supplied by mantle convection, Q _melt is heat loss due to the eruption of upwelling mantle melt at the surface, Q _rad is heat generated by radioactive decay in the mantle, Q _cmb is heat lost from the core across the core-mantle boundary (CMB), Q _man is the secular heat lost from the mantle, Q _tidal is heat generated in the mantle by tidal dissipation, and Q _L,man is latent heat released by the solidification of the mantle. Crustal heat sources have been excluded because they do not contribute to the mantle heat budget. Note that heat can be released from the mantle in two ways: via conduction through the upper mantle thermal boundary layer (Q _conv) and by melt eruption (Q _melt). Detailed expressions for heat flows and temperature profiles as functions of mantle and core properties are given in Appendix B.

The thermal evolution model assumes a mobile-lid style of mantle heat loss where the mantle thermal boundary layers maintain Rayleigh numbers that are critical for convection. In contrast, a stagnant-lid mantle parameterization would have a lower heat flow than a mobile lid at the same temperature (e.g., Solomatov and Moresi, 2000). However, a stagnant-lid mantle that erupts melt efficiently to the surface can lose heat as efficiently as a mobile-lid mantle with no melt heat loss (Moore and Webb, 2013; Driscoll and Bercovici, 2014). Io is an example of this style of mantle cooling (O'Reilly and Davies, 1981; Moore et al., 2007).

Similarly, the thermal evolution of the core is governed by the conservation of energy in the core, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}Q_{\rm cmb} = Q_{\rm core} + Q_{\rm icb} + Q_{\rm rad , core} \tag{11}\end{align*} \end{document}

where Q _core is core secular cooling, Q _rad,core is radiogenic heat production in the core, and heat released by the solidification of the inner core is Q _icb = \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\dot {M}$$ \end{document} _ic(L _Fe + E _G), where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\dot {M}_{ \rm ic}$$ \end{document} is the change in inner core mass M _ic, and L _Fe and E _G are the latent and gravitational energy released per unit mass at the inner-core boundary (ICB).

The internal thermal evolution equations are derived by using the secular cooling equation Q_i = −c_iM_i \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\dot {T}$$ \end{document} _i, where c is specific heat and i refers to either mantle or core, in Eqs. 10 and 11. Solving for \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\dot {T}$$ \end{document} _m and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\dot {T}$$ \end{document} _c gives the mantle and core thermal evolution equations, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\dot{T}_{\rm m} = ( Q_{\rm cmb} + Q_{\rm rad} + Q_{\rm tidal} + Q_{\rm L , man} - Q_{\rm conv} - Q_{\rm melt} ) / M_{\rm m} c_{\rm m} \tag{12}\end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\dot { T } _ { \rm c } = - \frac { ( Q_ { \rm cmb }- Q_{\rm rad , c } ) } { M_ { \rm c } c_ { \rm c } - A_ { \rm ic} \rho_ { \rm ic } \epsilon_ { \rm c } \frac { dR_ { \rm ic } } { dT_ { \rm cmb } } ( L_ { \rm Fe } + E_ { \rm G } ) } \tag { 13 } \end{align*} \end{document}

where the denominator of (13) is the sum of core specific heat and heat released by the inner core growth, A _ic is inner core surface area, ρ _ic is inner core density, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\epsilon_{ \rm c}$$ \end{document} is a constant that relates average core temperature to CMB temperature, dR _ic /dT _cmb is the rate of inner core growth as a function of CMB temperature, and L _Fe and E _G are the latent and gravitational energy released at the ICB per unit mass (Table 1). See Appendix B and Driscoll and Bercovici (2014) for more details.

2.3. Orbital evolution model

The orbital evolution of the planet's semimajor axis a and eccentricity e, assuming no dissipation in the primary body (the star), is (Goldreich and Soter, 1966; Ferraz-Mello et al., 2008; Jackson et al., 2009) \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\dot{e} = \frac {21} {2} Im ( k_2 ) \frac {M_*}{M_{\rm p}} \left(\frac{R_{\rm p}}{{a}}\right) ^5 ne \tag{14}\end{align*} \end{document}

and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\dot {a} = 2ea \dot{e} \tag{15}\end{align*} \end{document}

Mean motion can be replaced by using n² = \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$GM_*/ a^3$$ \end{document} , \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\dot { e } = \frac { 21 } { 2 } Im ( k_2 ) \frac { M_*^ { 3 / 2 } G^ { 1 / 2 } R_ { \rm p } { ^5 } } { M_ { \rm p } } \frac { e } { a^ { 13 / 2 } } \tag { 16 } \end{align*} \end{document}

The differential equations for thermal evolution (12, 13) and orbital evolution (16, 15) are solved simultaneously to compute coupled thermal-orbital evolutions.

3. Steady-State Solutions

Before exploring the full model, it is useful to highlight the influence of tidal heating on the thermal evolution in a steady-state sense by comparing heat flows as functions of mantle temperature. Figure 3 shows the tidal heat flow (a) and orbital circularization time (b), \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}t_{\rm circ} = e / \dot{e} \tag{17}\end{align*} \end{document}

as a function of mantle temperature for a range of initial orbital distances and eccentricity of e = 0.1. Figure 3a shows that a peak in dissipation occurs when the mantle is in a partially liquid viscoelastic state (T _m ≈ 1800 K), where initial tidal perturbations behave elastically and the long timescale relaxation is viscous. Dissipation is lower in a colder mantle where the response is closer to purely elastic. Dissipation is also lower in a hotter, mostly liquid mantle where the material behaves viscously, with little resistance to the external forcing.

FIG. 3.

Tidal dissipation properties as a function of mantle temperature T _m for M _* = 0.1M _sun and e = 0.1. (a) Comparison of tidal heat flow from the Maxwell model (curves) with the constant \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\cal {Q}$$ \end{document} model ( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\cal {Q}$$ \end{document} = 100, k ₂ = 0.3) at four orbital distances (see legend). Also shown in (a) is the mantle surface heat flow Q _surf as a function of temperature (solid black) and constant runaway greenhouse threshold (dashed). (b) Timescales for orbital circularization using the Maxwell model [same colors as in (a)] and mantle cooling (solid black).

Also plotted in Fig. 3a is the convective mantle cooling curve from (24), which reflects the preferred cooling rate of the interior. Conceptually, a tidal steady state is achieved as the planet cools down from an initially hot state (T _m >2000 K) until the convective cooling curve intersects the tidal dissipation (heat source) rate. This intersection implies that heat loss is in balance with tidal heating so that the interior stops cooling. The steady state occurs around 1850–1950 K over the range of orbital distances in Fig. 3a. The steady state is maintained until the orbit begins to circularize, which drops the dissipation curve and intersection point to lower temperatures. Circularization continues slowly until the dissipation rate falls below the surface cooling rate, at which point the planet resumes cooling normally with tidal heat playing a minor role in the heat budget.

The time required to circularize, shown in Fig. 3b, is inversely proportional to dissipation rate through (16 and 17) so that a mantle in a viscoelastic state dissipates orbital energy efficiently and circularizes quickly. At the inner edge of the HZ (a = 0.02 AU), circularization occurs in less than ∼1 Gyr, while on the outer edge circularization requires billions of years or may not occur at all. Also shown in Fig. 3b is the mantle cooling time t _T = M _m c _m(T _m(0) − T _m)/Q _surf, which is the time required for the mantle to cool from T _m(0) = 2500 K to T _m. This shows that it takes the mantle ∼1 Gyr to adjust to a change in the tidal heat source. The cooling time is typically longer than the circularization time in the HZ (Fig. 3b), implying that tidal heating can evolve faster than the thermal response of the mantle.

4. Model Results: Evolving Orbits

This section presents full thermal-orbital evolutions with a single Earth-mass planet in orbit around a 0.1 solar mass M star. We first focus on planets orbiting a 0.1 solar mass star, where the HZ is very close to the star, in order to examine an extreme tidal environment for a lone planet. Section 7 investigates a range of stellar masses. In this section the orbit of the planet is free to evolve according to Eqs. 16 and 15. Later in Section 5 we explore thermal evolutions with fixed orbits. The models all have A_μ = 2 × 10⁵ J mol⁻¹, T _m(0) = 2400 K, and T _c(0) = 6000 K. The results are independent of initial mantle and core temperatures up to approximately ±500 K.

4.1. Influence of initial orbital distance

First, we investigate the evolution of three models that start with e(0) = 0.5 at three orbital distances: (1) inside the inner edge of the HZ at a = 0.01 AU, (2) within the HZ at a = 0.02 AU, and (3) the outer edge of the HZ at a = 0.05 AU.

Figure 4 compares the tidal heat flow and eccentricity of these three models as a function of T _m. The inner case with a(0) = 0.01 AU begins with a rather high initial tidal heat flow (Q _tidal∼0.1 TW), considering the mantle is mostly molten at this time. Being so close to the star, the planet has a fast circularization time (Fig. 3), so eccentricity decreases rapidly and tidal dissipation effectively ends within 1 Myr (also see Fig. 5c). This implies that circularization occurs during the magma ocean stage.

FIG. 4.

Thermal-orbital evolution for three models with initial orbits of e(0) = 0.5 and a(0) = 0.01 (red), 0.02 (black), and 0.05 (blue). The temperature that corresponds to 50% melt fraction is denoted by the dashed line.

FIG. 5.

Time evolution of models with initial orbits of e = 0.5 and a = 0.01 (red), 0.02 (black), and 0.05 (blue). (a) Temperature in the mantle (solid) and core (dashed). (b) Heat flow at the surface Q _surf (solid), tidal dissipation Q _tidal (dash-dot), mantle radiogenic heating (dotted), mantle latent heat (dash-dot-dot), and core heat flow Q _cmb (dashed). The runaway greenhouse heat flow threshold (1.53 × 10⁵ TW) is labeled as a solid gray line. (c) Orbital distance. (d) Eccentricity. (e) Magnetic moment of core dynamo (solid) and inner core radius (dashed). Inner core radius axis has been scaled by core radius so the top corresponds to a completely solid core. For reference, Earth's present-day magnetic moment is about 80 ZAm². (f) Melt mass flux to the surface. Melt eruption fluxes for present-day mid-ocean ridges (10¹³ kg yr⁻¹) and the Siberian traps (10¹⁵ kg yr⁻¹) shown for reference (gray dashed).

The middle case with a(0) = 0.02 AU begins with a lower tidal heat flow because it orbits farther from the star. As the mantle cools and solidifies, tidal dissipation evolves through the viscoelastic peak at T _m≈1800 K where the mantle is ∼50% molten and a peak of Q _tidal∼100 TW occurs. This increase in dissipation drives a rapid circularization (Fig. 4b), which then decreases the dissipation as the mantle cools further.

The outer case at a(0) = 0.05 AU experiences the lowest initial tidal heat flow of Q _tidal(0)∼10⁻⁵ TW due to it being farthest from the star. Dissipation remains low, and the orbit remains eccentric until the mantle cools to T _m∼1800 K, at which point dissipation increases rapidly. Tidal heat flow peaks around Q _tidal∼100 TW and T _m∼1750 K before decreasing due to decreasing eccentricity. The peak in dissipation occurs at a slightly lower temperature than the middle case because the eccentricity remains higher longer due to slower circularization. In fact, the circularization time is slow enough that after 10 Gyr the model still retains a finite eccentricity of e∼0.01 and a tidal heat flow of Q _tidal∼ 10 TW. This shows that tidal dissipation at the outer edge can linger longer due to slower circularization times.

A detailed comparison of these three models over time is shown in Fig. 5. Relatively small differences in their temperature histories (Fig. 5a) are driven by small differences in mantle and core cooling rates (Fig. 5b). Circularization of the inner model occurs in the first million years and by 100 Myr for the middle model, while the outer model retains a small eccentricity of e = 2 × 10⁻³ after 10 Gyr (Fig. 5d). These circularization times are reflected in the tidal heat flow peaks (Fig. 5b), which occur around 0.1 Myr for the inner case, 10 Myr for the middle case, and 1 Gyr for the outer case. Inward migration by 10–20% also accompanies this circularization (Fig. 5c).

The thermal evolutions are mainly controlled by secular cooling and radiogenic decay, with tidal dissipation as a temporary energy source. Mantle heating due to latent heat released during the solidification of the mantle is of order 10⁴ TW until ∼1 Myr, then drops below ∼1 TW once the mantle is mostly solid around 0.1 Gyr. This decrease in latent heat causes the mantle heat flow to drop rapidly between 1 and 10 Myr. Mantle solidification and the drop in mantle heat flow occur slightly later for the inner model for two reasons: (1) the surface heat flow is lower than the other models in the first million years because the surface is hot, decreasing the upper mantle temperature jump; (2) tidal heating is initially moderate (∼0.1 TW) despite the mantle being mostly molten due to proximity to the star.

These mobile-lid Earth-like models have a strong temperature feedback, or thermostat effect, such that if mantle temperature increases (e.g., due to tidal dissipation), the viscosity decreases rapidly, and the boundary layers thin out, resulting in an increase in the boundary heat flows. Consequently, increases in internal heat sources are accommodated by increases in heat flows such that the mantle and core cool monotonically. One minor exception is the brief heating of the core at 1 Gyr due to early radioactive decay in the core. In contrast, a stagnant-lid parameterization with a weaker heat flow–temperature feedback would force the mantle to maintain higher temperatures and thus accommodate the same cooling rates (e.g., Solomatov and Moresi, 2000; Driscoll and Bercovici, 2014). Therefore, we expect mobile-lid planets to cool faster, dissipate tidal energy more efficiently, and circularize faster than stagnant-lid planets. For stagnant-lid planets that are strongly tidally heated, melting rather than conduction removes heat from the interior, as Io demonstrates today.

Core cooling rates are similar at these three orbital distances, which results in similar magnetic moment histories and inner core nucleation times (Fig. 5e). Inner core nucleation induces a kink in the core compositional buoyancy flux and magnetic moment around 4 Gyr, similar to predictions for Earth (Driscoll and Bercovici, 2014). Surface melt eruption rate is determined by the mantle cooling rate through the upper mantle geothermal gradient, so that the eruption rates at these three orbital distances are similar and follow the mantle heat flow history (Fig. 5f). After 6 Gyr, the middle and outer planets' mantles are completely solid so that melt eruption ends, while melt eruption continues longer for the inner case due to a slightly hotter mantle.

In summary, Earth-like planets near the inner edge of the HZ around 0.1 M _* stars circularize rapidly (within a few million years), allowing internal cooling and core dynamo action to proceed similar to the Earth model. On the outer edge of the HZ, orbital circularization is slower, which leads to a prolonged period of tidal dissipation that is accentuated by the cooling of the mantle through a viscoelastic state after ∼1 Gyr. Despite these differences in the tidal evolution, the magnetic and magmatic evolutions of these mobile-lid planets are similar. These three cases with high initial eccentricities of e(0) = 0.5 demonstrate the potential for a strong coupling between orbital and thermal evolution by tidal dissipation.

4.2. Summary evolution contours

In this section, we compare the final states (after 10 Gyr) of orbital-thermal evolution for 132 models that span a range of initial orbital distances of 0.01–0.10 AU and initial eccentricities of 0–0.5. The results are displayed as contours in initial orbital a-e space (Figs. 7 –13). Here, we consider models whose orbits evolve (left panels of Figs. 7 –13), while Section 5 below considers models whose orbits are fixed (right panels of Figs. 7 –13).

Figure 6 shows the fractional change in orbital distance (a) and eccentricity (b). After 10 Gyr, most models have circularized within the HZ due to tidal dissipation (Fig. 6b). The iso-contour lines in Fig. 6b are nearly vertical because eccentricity evolution is proportional to e/a ^13/2 in (16). In other words, orbital circularization is a stronger function of orbital distance than eccentricity. Circularization also causes the orbit to migrate inward (Fig. 6a), although this results in a maximum inward migration of only 22% of the initial distance. The evolution of orbital distance, which produces mainly horizontal iso-contours (Fig. 6a), is controlled by initial eccentricity because migration is proportional to e ²/a ^11/2 by (15); hence, migration \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$( \dot {a} )$$ \end{document} is a stronger function of eccentricity than circularization \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$( \dot {e} )$$ \end{document} .

FIG. 6.

Contour of orbital evolution after 10 Gyr for a range of initial orbital distances and eccentricities. (a) Change in orbital distance: (a − a ₀)/a ₀. (b) Change in eccentricity: (e − e ₀)/e ₀. Orbits are free to evolve in both panels. The HZ in denoted by vertical dashed white lines.

Figure 7a contours tidal heat flow for these models. We identify the tidal heat flow boundaries defined by Barnes and Heller (2013) as an Earth Twin for Q _tidal < 20 TW, a Tidal Earth for 20 < Q _tidal < 1020 TW, and a Super-Io for Q _tidal >1020 TW. Models that circularize by 10 Gyr have zero tidal heating. At the outer edge of the HZ, circularization is still occurring at a rate that is proportional to the change in e in Fig. 6b. Beyond a∼0.07 AU, tidal dissipation is too weak to result in any significant circularization. Therefore, there are gradients in Q _tidal on both sides of this boundary at a∼0.07. There is also a decrease in Q _tidal with initial e because models with low initial e circularize earlier. The combination of these effects results in a peak in tidal dissipation around a∼0.07 AU and e∼0.5.

FIG. 7.

Contour of (log) tidal heat flow after 10 Gyr for a range of initial orbital distances and eccentricities. (a) Orbit evolves. (b) Orbit is fixed. The tidal heat flow boundaries defined by Barnes and Heller (2013) are shown for Earth Twins Q _tidal < 20 TW, Tidal Earths 20 < Q _tidal < 1020 TW, and Super-Ios for Q _tidal > 1020 TW.

This peak in tidal heat flow causes a slight increase in heat flow (Fig. 8a) and mantle temperature (Fig. 10a). This is an example of the steady-state behavior of Fig. 3 where higher tidal heat flows intersect the convective heat flow curve at higher mantle temperatures. Hotter mantle temperatures also cause the lower mantle to have lower viscosity, thinner boundary layers, and increased core heat flows (Fig. 9a). A second maximum in core heat flow occurs at the innermost orbits due to the high surface temperature insulating the mantle, which keeps mantle temperature high (Fig. 10a) and thins the lower mantle thermal boundary layer. Core temperature is low where core heat flow is high (Fig. 11a) due to secular cooling of the core.

FIG. 8.

Contour of surface heat flow after 10 Gyr for a range of initial orbital distances and eccentricities. (a) Orbit evolves. (b) Orbit is fixed. White contour line shown at Earth's present-day surface heat flow ( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$Q^*_{ \rm surf}$$ \end{document} = 40 TW).

FIG. 9.

Contour of core heat flow after 10 Gyr for a range of initial orbital distances and eccentricities. (a) Orbit evolves. (b) Orbit is fixed.

FIG. 10.

Contour of mantle temperature after 10 Gyr for a range of initial orbital distances and eccentricities. (a) Orbit evolves. (b) Orbit is fixed.

FIG. 11.

Contour of core temperature after 10 Gyr for a range of initial orbital distances and eccentricities. Line contours show solid core fraction R _ic/R _c as a percentage (i.e., 100% corresponds to a completely solid core). (a) Orbit evolves. (b) Orbit is fixed.

After 10 Gyr of significant core cooling, all models have a large solid inner core. The size of the inner core (also contoured in Fig. 11) is proportional to Q _cmb (Fig. 9a) and is between 80% and 100% of the core radius (also see Fig. 5e). Where the core is entirely solid, no dynamo action is possible (upper left corner of Fig. 12a), and where the core is mostly solid, the magnetic moment is weak due to the small size of the dynamo region (upper right corner of Fig. 12a).

FIG. 12.

Contour of magnetic moment after 10 Gyr for a range of initial orbital distances and eccentricities. (a) Orbit evolves. (b) Orbit is fixed. For reference, Earth's present-day magnetic moment is about 80 ZAm².

The eruption of melt to the surface (Fig. 13a) is controlled by the upper mantle geothermal gradient and thus proportional to mantle heat flow with a peak around 0.07 AU. A secondary peak in melt mass flux at close-in orbits (upper left corner of Fig. 13a) is caused by a slightly higher mantle temperature associated with a hotter, insulating surface (Fig. 10a).

FIG. 13.

Contour of surface melt mass flux after 10 Gyr for a range of initial orbital distances and eccentricities. (a) Orbit evolves. (b) Orbit is fixed. White line contour denotes Earth's approximate present-day mid-ocean-ridge melt flux (10¹³ kg yr⁻¹). Note color scales in (a) and (b) are different.

5. Model Results: Fixed Orbits

In this section, we consider planets whose orbits are fixed \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$( \dot {e} = \dot {a} = 0 )$$ \end{document} . This includes eccentric orbits, which could be fixed, for example, through interactions with a planetary companion (Van Laerhoven et al., 2014). Figures 7b –13b show contours in orbital space, similar to those discussed above except with fixed orbits.

5.1. Example of tidal steady state

Figure 14 shows the time evolution of a specific case with a fixed orbit of a = 0.02 AU and e = 0.5 that reaches a tidal steady state, an example of the scenario described in Section 3. Tidal heating initially starts low (∼10⁻³ TW) before increasing as the mantle cools for the first 10 Myr, until the mantle reaches a steady-state temperature of T _m ≈ 1800 K. In this steady state, heat loss is balanced by internal sources so that mantle cooling becomes insignificant. The steady-state surface heat flow (Q _surf ≈ 1000 TW) corresponds to a surface heat flux similar to that of Io (q _surf ≈ 2 W m⁻²), implying that this planet might be better characterized as a super-Io than Earth-like (e.g., Barnes et al., 2010). The tidal steady state still allows the core to cool slowly because a significant temperature difference between the mantle and core persists. Core cooling drives a core dynamo for all 10 Gyr, although the magnetic moment rapidly declines as the core is nearly entirely solid by 10 Gyr (Fig. 14c).

FIG. 14.

Time evolution of a model with fixed orbit of e = 0.5 and a = 0.02 AU. (a) Temperature in the mantle (solid) and core (dashed). (b) Heat flow at the top of the mantle Q _conv (solid), tidal Q _tidal (dash-dot), mantle radiogenic heating (dotted), and core heat flow Q _cmb (dashed). (c) Magnetic moment of core dynamo (solid) and inner core radius (dashed). Inner core radius axis goes from zero to total core radius. For reference, Earth's present-day magnetic moment is about 80 ZAm². (d) Melt mass flux to the surface. Melt eruption fluxes for present-day mid-ocean ridges (10¹³ kg yr⁻¹) and the Siberian traps (10¹⁵ kg yr⁻¹) shown for reference (gray dashed).

5.2. Summary contours

These fixed-orbit models, in contrast with the evolving models in Section 4, have tidal heat flows that are mainly determined by the orbital state rather than the cooling history. Specifically, Q _tidal increases with e and decreases with a, producing a maximum in the upper left corner of Fig. 7b. Mantle heat flow (Fig. 8b) and temperature (Fig. 10b) increase with tidal heating due to the positive feedback between mantle temperature and surface heat flow.

Core heat flow peaks in models at moderate orbital distances where tidal heat flow is similar in magnitude to the sum of all other mantle heat sources, that is, Q _tidal∼Q _cmb +Q _man + Q _rad (Fig. 9b). This peak can be understood by considering how Q _cmb behaves at the two tidal extremes: (1) where tidal heating is strong (upper left corner of Fig. 7b) the mantle is forced into a hot steady state so that surface heat flow can accommodate all heat sources, which thins the lower mantle thermal boundary layer and allows a moderate core heat flow of Q _cmb∼10 TW; (2) where tidal heating is weak (lower right corner of Fig. 7b) the mantle and core are free to cool similar to the Earth model, so that Q _cmb decreases monotonically over time. In between these limits, tidal dissipation heats the mantle slightly, increasing the surface heat flow, but does not dominate the heat budget, which allows the interior to cool. Note that even when mantle temperature is high (∼2000 K) it is still significantly colder than the core (∼3800 K) so that the mantle and core are not in thermal equilibrium and the core is forced to cool. The effect is to produce a region where a modest amount of tidal dissipation actually promotes core cooling, similar to the peak in the evolving models (Fig. 9a). We refer to this 30% increase in Q _cmb as the super-cooling of the core.

The influence of this peak in core cooling rate on the dynamics of the core is dramatic. Intuitively, core temperature is lowest where core cooling is high (center of Fig. 11b), and the coldest models with T _cmb ≈ 3850 K have completely solid cores (i.e., R _ic = R _core). A fully solid core prevents fluid motion and therefore dynamo action. These models lose their dynamos at ∼9.8 Gyr, so they have dynamos a vast majority of the time. Note that our magnetic scaling law likely provides an upper limit on the dynamo lifetime because the scaling law was derived from thick shell dynamos and does not account for stratified layers.

This prediction implies that there is a dip (and possibly a gap) in magnetic field strength for tidally heated and orbitally fixed planets over most of the HZ (Fig. 12b). We note that our core liquidus does not include light element depression (e.g., Hirose et al., 2013), which would tend to slow inner core nucleation and allow some of these models to maintain a liquid region slightly longer. This result emphasizes the difference between core cooling and dynamo action; cooling is ongoing (at least until thermal equilibration), whereas dynamo lifetime, which relies on convection in the liquid, is limited by the solidification time of the core (see also Gaidos et al., 2010; Tachinami et al., 2011). In other words, rapid core cooling is helpful for temporarily driving dynamo action but shortens the lifetime of the dynamo.

The eruption of mantle melt to the surface follows surface heat flow (Fig. 13b). The extreme mass fluxes of 10¹⁶ kg yr⁻¹ correspond to a global basalt layer resurfacing rate of 7 m kyr⁻¹. For reference, the Siberian traps, one of the largest igneous provinces on Earth and thought to be responsible for the Permian mass extinction event, is estimated to have produced a basalt layer at a rate of 1 m kyr⁻¹ over the area of the traps (Reichow et al., 2002). Therefore, continuous eruption rates of ∼10¹⁶ kg yr⁻¹ are likely to prevent such planets from being habitable.

We also compute the same range of models but only fixing eccentricity, allowing a to evolve. This might occur if a neighboring planet forces the eccentricity but allows inward migration. In these cases, we find that all models with initial orbits of a < 0.05 AU (or a < 0.02 AU) and e > 0.2 (or e > 0.1) migrate into the central star within 10 Gyr, and most by 5 Gyr.

6. Internally Driven Runaway Greenhouse

As described by Barnes et al. (2013), if interior heat flux exceeds the limit at which energy can be radiated from the top of the atmosphere, then runaway heating of the surface occurs, which evaporates the ocean and leads to rapid water loss (Goldblatt and Watson, 2012). Figure 15 shows the time spent in an internally driven runaway greenhouse, defined as the period of time when the surface heat flux exceeds the threshold q _runaway = 300 W m⁻².

FIG. 15.

Contour of time spent in an internally driven runaway greenhouse, defined as when surface heat flow exceeds the threshold for a runaway greenhouse (300 W m⁻², or 1.53 × 10¹⁷ W), for a range of initial orbital distances and eccentricities. (a) Orbit evolves. (b) Orbit is fixed.

For both evolving (Fig. 15a) and fixed (Fig. 15b) orbital models, the runaway greenhouse period is shorter at closer orbital distances, almost independent of eccentricity. This implies that tides, which depend strongly on eccentricity, play a minor role in the length of the runaway greenhouse state. The runaway greenhouse state is shorter for close-in planets because they have higher effective surface temperatures closer to the star, which insulates the mantle and decreases the initial surface heat flow. With lower initial surface heat flows, these inner planets drop below the runaway heat flow threshold earlier (Fig. 5b). A second trend in Fig. 15a toward even shorter times spent in a runaway greenhouse is found for the innermost, high-eccentricity planets. This drop in surface heat flow at around 50–100 kyr occurs during the circularization of the inner planets' orbits, when the tidal heat flow rapidly declines (Fig. 5b). Circularization causes a small dip in the surface heat flow as the interior temperatures and heat flows adjust to the smaller internal (tidal) heat source. This adjustment to lower heat flows, although seemingly minor, actually shortens the time spent above the runaway threshold (Fig. 5b).

Interestingly, when the heat flow is high enough to drive a runaway greenhouse during the first few hundred million years, the mantle is so hot that tidal dissipation is inefficient. Typically tidal dissipation is not a major heat source until the mantle solidifies and cools down to ∼1800 K, which occurs after the runaway greenhouse and magma ocean phases.

In summary, we find that mobile-lid Earth-like planets typically spend several hundred thousand years in an internally driven atmospheric runaway greenhouse state and that tidal dissipation in the mantle at this time plays a minor role. The runaway greenhouse timescale (∼100 kyr) is shorter than the typical magma ocean solidification time (∼10 Myr), a period when the surface is likely uninhabitable anyway. These models assume mobile-lid cooling at all times; however, Foley et al. (2012) proposed that a runaway greenhouse could induce a transition from mobile to stagnant lid, which would also slow internal cooling and be detrimental to habitability. In Section 7, we explore more generally how tides may affect habitability by computing the length of time spent in a tidally dominated state for a range of stellar masses.

7. Influence of Stellar Mass

The above calculations assumed a stellar mass of 0.1 M _sun. In this section, we explore the influence of stellar mass, in the range 0.1–0.6 M _sun, at the inner edge of the HZ. Similar to the contours in Section 4.2, we compute a grid of models with a range of initial eccentricities of 0–0.5 and allow the orbit to evolve in time. The initial orbital distance is set just outside the inner edge of the radiative HZ, which is derived from the stellar mass by the parametric equations of Kopparapu et al. (2014), so that the planet remains in the HZ after 10 Gyr of orbital migration.

Figure 16 summarizes the results of these models in terms of two timescales: (a) the time spent in a tidally dominated state, defined as when the tidal heat flow Q _tidal is 50% or more of the total surface heat flow Q _surf; (b) the time to reach Earth's present-day surface heat flow of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$Q^*_{ \rm surf} = 40 \ { \rm TW}$$ \end{document} .

FIG. 16.

Contour of (a) time spent in a tidally dominated state (i.e., Q _tidal/Q _total ≥ 0.5) and (b) time to reach Earth's present-day surface heat flow (Q _surf = 40 TW). In (b) a white contour line is shown at 4.5 Gyr.

The islandlike shapes of these time contours can be explained by a combination of three physical effects. First, planets with initially low eccentricity (e(0) < 0.1) experience weak tides and spend little time, if any, in the tidally dominated regime. At higher eccentricity, tides become stronger, so that more eccentric planets are tidally dominated longer (Fig. 16a). Second, as stellar mass increases, the HZ moves to larger orbital distances, and the tidal dissipation decreases because tidal dissipation in (2) is a stronger function of orbital distance (∝ a ^−15/2) than stellar mass \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$( \propto M_*^{ + 5 / 2} )$$ \end{document} . The net result is a decrease in tidal dissipation within the HZ for increasing stellar mass, and a shorter time spent in the tidally dominated state (Fig. 16a). This effect produces contour boundaries with positive slope in Fig. 16. Third, models with high initial eccentricity (e(0) > 0.2) and close-in initial orbits around low-mass stars (M _star < 0.12) experience extreme early tides that drive rapid orbital circularization. This leads to short times spent in the tidally dominated state.

Figure 16b, similar to Fig. 16a, shows that eccentric planets on the inner edge around 0.15–0.4 M _sun stars maintain surface heat flows in excess of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$Q^*_{ \rm surf}$$ \end{document} for 10 Gyr due to strong tidal dissipation. Interestingly, Fig. 16b shows that planets that experience only a temporary period of tidal heating actually cool to an Earth-like heat flow before 4.5 Gyr. These planets cool faster than Earth because their thermal adjustment timescale is longer than their circularization (or tidal heating) timescale, so they are still adjusting to the new heat balance with a lower tidal heat source. In other words, the surface heat flow that was increased during the tidal heating phase is still slightly larger than it would have been with no tidal heating. This super-cooling effect was also discussed in Section 4.2.

In summary, tides are more influential around low-mass stars. For example, planets around 0.2 M _sun stars with eccentricity of 0.4 experience a tidal runaway greenhouse for 1 Gyr and would be tidally dominated for 10 Gyr. These timescales would increase if the orbits were fixed, for example, by perturbations by a secondary planetary companion. We find a threshold at a stellar mass of 0.45 M _sun, above which the HZ is not tidally dominated. These stars would be favorable targets in the search for geologically habitable Earth-like planets, as they are not overwhelmed by strong tides.

8. Discussion

In summary, we have investigated the influence of tidal dissipation on the thermal-orbital evolution of Earth-like planets around M stars with masses 0.1–0.6 M _sun. A thermal-orbital steady state is illustrated where, under certain conditions, heat from tidal dissipation is balanced by surface heat flow. We find that mantle temperatures in this balance are hotter for planets with shorter orbital distances and larger eccentricities. Orbital energy dissipated as tidal heat in the interior drives both inward migration and circularization, with a circularization time that is inversely proportional to the dissipation rate. The cooling of an eccentric planet in the HZ leads to a peak in the dissipation rate as the mantle passes through a viscoelastic rheology state. Planets around 0.1 solar mass stars with initial orbits of a < 0.07 AU circularize before 10 Gyr, independent of initial eccentricity. Once circular, these planets cool monotonically and maintain dynamos similar to that of Earth. Generally, we find that tidal dissipation plays a minor role on the dynamo history if the orbit is free to evolve in time.

When the orbit is fixed, the planet cools until a tidal steady-state balance between tidal dissipation and surface cooling is reached. In the HZ, this steady state can produce a super-cooling of the core when tidal heating is strong enough to heat the mantle and decrease its viscosity and low enough not to dominate the surface heat flow. This rapid cooling leads to complete core solidification, prohibiting dynamo action for most models in the HZ with e > 0.05 by 10 Gyr. In addition to weak magnetic fields, massive melt eruption rates in the HZ may render these fixed-orbit planets uninhabitable.

Commonly, the term “habitability” refers to the influx of radiation necessary to maintain surface liquid water. However, the full habitability of a planet must involve the dynamics of the interior and its interaction with the surface environment. We find that tidal heating of a planetary mantle can influence surface habitability in several important ways:

(1) Prolonged magma ocean stage. Close-in planets with a high eccentricity (e ≳ 0.1) will experience extreme tidal heating rates of ∼1000 TW and tidal steady-state mantle temperatures of ∼2000 K, implying mostly molten mantles. These super-tidal planets are uninhabitable, as the surface itself is likely molten or close to the silicate solidus.

(2) Extreme volcanic eruption rates. Tidal heating, in addition to increasing surface heat flow, can produce extreme surface melt production rates. Even if only a fraction (∼20%) of this melt erupts to the surface, it can easily produce a 100-fold increase over the present-day mid-ocean-ridge eruption rate (∼10¹³ kg yr⁻¹). These extreme eruption rates can lead to rapid global resurfacing and degassing that render the surface environment a violent and potentially toxic place for life. Volcanically dominated atmospheres could be significantly different from modern-day Earth's and are potentially detectable with future space- and ground-based telescopes (Misra et al., 2015).

(3) Lack of magnetic field. Planetary magnetic fields are often invoked as shields necessary to maintain life. Magnetic fields can protect the atmosphere from stellar wind erosion (Driscoll and Bercovici, 2013) and the surface from harmful radiation (Griessmeier et al., 2005; Dartnell, 2011). Super-cooling of the core, which can solidify the entire core and kill the dynamo, occurs in the HZ after ∼9 Gyr with a fixed orbit. Alternatively, a tidally heated stagnant-lid planet can maintain hotter mantle temperatures and lower core cooling rates, weakening the core-generated magnetic field. Even before losing the dynamo entirely, these planets may have magnetic fields that are too weak to hold the stellar wind above the atmosphere or surface. In either case, the lack of a strong magnetic shield will be detrimental to life.

(4) Tidally driven runaway greenhouse. In Section 6, we show that driving a runaway greenhouse by tidal heating in the rocky interior alone is difficult. To achieve the runaway threshold heat flux at the surface, either the planet would have to be forced into a highly eccentric orbit after the mantle has cooled down to ∼1800 K, or the dissipative material properties would have to be different. For example, if the mantle were composed of a lower-viscosity material, then the maximum Maxwell tidal power could increase to 10⁵ TW (Fig. 2a). A significant amount of tidal energy can also be dissipated in the liquid portions of the planet (Tyler, 2014), which is beyond the scope of this study.

With growing interest in the habitability of Earth-like exoplanets, the development of geophysical evolution models will be necessary to predict whether these planets have all the components that are conducive for life. This paper focused on a single Earth-mass planet, but the mathematical equations can be developed to model the evolution of other rocky planet/star mass ratios, including large rocky satellites around giant planets. However, significant uncertainties make the application to super-Earths particularly challenging. The fundamental physical mechanisms underpinning plate tectonics, both in terms of its generation and maintenance over time, are not fully understood, which makes extrapolation to larger planets questionable. Perhaps most importantly, material properties, such as viscosity, melting point, solubility, and conductivity, are poorly constrained at pressures and temperatures more extreme than Earth's lower mantle and core. This uncertainty prevails in our own solar system where the divergent evolution of Earth and Venus from similar initial conditions to dramatically different present-day states remains elusive.

Future thermal-orbital modeling improvements should include coupling the evolution of the interior to the surface through volatile cycling and atmosphere stability. Advancing the orbital model to include gravitational interactions with additional planetary companions would allow for tidal resonances, variable rotation rates, and other time-dependent orbital forcings. In addition to the eccentricity tide explored here, an obliquity tide could also be important. Further improvements could include dissipation in oceans or internal liquid layers, variable internal composition, internal structures, radiogenic heating rates, core light element depression, continental crust formation, and eventually a direct coupling of first-principles numerical simulations.

Appendix A. Tidal Dissipation Model

This section demonstrates the dependence of the tidal dissipation model and material properties on mantle temperature. Figure A1 shows several parameters related to the tidal dissipation rate as a function of mantle temperature for the nominal shear modulus activation energy of A_μ = 2 × 10⁵ J mol⁻¹. The Maxwell model that uses the full form of −Im(k ₂) in (4) differs from the common approximation of −Im(k ₂) ≈ k ₂/ \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\cal {Q}$$ \end{document} for mantles hotter than the present-day (T _m >1630 K) (Fig. A1c). The difference between the Maxwell model and this approximation corresponds to about 10 orders of magnitude larger tidal heat flow at high temperature (Fig. A1d). The approximation is invalid at high temperature because it does not account for the drop in tidal dissipation expected in a liquid, since the approximation relies on \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\cal {Q} \propto \eta / \mu$$ \end{document} , which is constant, whereas the Maxwell model predicts a sharp drop in tidal dissipation with viscosity when μ/β << 2/10.

FIG. A1.

Tidal dissipation parameters as a function of mantle temperature T _m for the nominal shear modulus activation energy of A_μ = 2 × 10⁵ J mol⁻¹. Circles denote the chosen calibration points for the present-day mantle.

Appendix B. Thermal History Model

B.1. Geotherm

The mantle temperature profile is assumed to be adiabatic everywhere except in the thermal boundary layers where it is conductive. The adiabatic temperature profile in the well-mixed region of the mantle is approximated to be linear in radius, which is a good approximation considering that mantle thickness D = 2891 km is much less than the adiabatic scale height H = c_p /αg ≈ 12,650 km, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} T_{\rm ad} = T_{\rm{UM}} + \gamma_{\rm ad} (R - r-\delta_{\rm UM})\tag{18}\end{align*} \end{document}

where the adiabatic gradient is γ _ad ≈ 0.5 K km. In the thermal boundary layers, the conductive temperature solutions, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \Delta T_{\rm{UM}} {\rm erf} \left[\frac{R-r} \delta_{\rm{UM}}\right] + T_{\rm s} \ , \ { \rm Upper \ mantle} \tag{19} \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\Delta T_ { \rm LM } { \rm erf } \left[ \frac { R_ { \rm c } - r } { \delta_ { \rm LM } } \right] + T_ { \rm cmb } \ , \ { \rm Lower \ mantle } \tag { 20 } \end{align*} \end{document}

replace the adiabat. Thermal boundary layer temperature jumps are ΔT _UM = T _UM − T _g and ΔT _LM = T _cmb − T _LM, and thermal boundary layer depth is δ. Figure B1 shows an example whole-planet geotherm T(r) at four times in the evolution. Surface temperature T _g is assumed to be equal to the equilibrium temperature, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}T_ { \rm eq } = \left( \frac { L_* } { 16 \pi \sigma a^2 } \right) ^ { 1 / 4 } \tag { 21 } \end{align*} \end{document}

where L _* is stellar luminosity and σ is the Stefan-Boltzmann constant.

The core temperature profile is assumed to be adiabatic throughout the entire core; that is, the thermal boundary layers within the core are ignored. This is a good approximation because the low viscosity and high thermal conductivity of liquid iron produce very small thermal boundary layers that are insignificant on the scale of the whole planet. The core adiabatic profile is approximated by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}T_ { \rm c } ( r ) = T_ { \rm cmb } \ { \rm exp } \left( \frac { R_ { \rm c } ^2 - r^2 } { D_N^2 } \right) \tag { 22 } \end{align*} \end{document}

where D_N ≈ 6340 km is an adiabatic length scale (Labrosse et al., 2001). The iron solidus is approximated by Lindemann's Law, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}T_ { \rm Fe } = T_{\rm { Fe , 0 }} \ { \rm exp } \left[ - 2 \left( 1 - \frac { 1 } { 3 \gamma_ { \rm c } } \right) \frac { r^2 } { D_ { \rm Fe } ^2 } \right] \tag { 23 } \end{align*} \end{document}

where T _{Fe, 0} = 5600 K, γ _c is the core Gruneisen parameter, and D _Fe = 7000 km is a constant length scale (Labrosse et al., 2001). This simple treatment of the core solidus does not account for volatile depression of the solidus, which has been demonstrated experimentally (Hirose et al., 2013), and would act to slow inner core growth. Inner core radius can then be solved for by finding the intersection of (22) and (23). For details, see Driscoll and Bercovici (2014).

B.2. Mantle and core heat flows

In this section, we define the remaining heat flows that appear in the mantle (10) and core (11) energy balance.

The convective cooling of the mantle Q _conv is proportional to the temperature gradient in the upper mantle thermal boundary layer, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}Q_ { \rm conv } = Ak_ { \rm UM } \frac { \Delta T_ { \rm UM } } { \delta_ { \rm UM } } \tag { 24 } \end{align*} \end{document}

where A is surface area and k _UM is upper mantle thermal conductivity. Q _conv is written in terms of T _m and the thermal boundary layer thickness δ _UM by requiring that the Rayleigh number of the boundary layer Ra _UM be equal to the critical Rayleigh number for thermal convection Ra _c ≈ 660 (Howard, 1966; Solomatov, 1995; Sotin and Labrosse, 1999; Driscoll and Bercovici, 2014). This constraint gives \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}Q_ { \rm conv } = Ak_ { \rm UM } \left( \frac { \alpha g } { Ra_ { \rm c } \kappa } \right) ^ { \beta } ( \epsilon_ { \rm UM } \Delta T_ { \rm m } ) ^ { \beta + 1 } ( \nu_ { \rm UM } ) ^ { - \beta } \tag { 25 } \end{align*} \end{document}

where the thermal boundary layer temperature jump ΔT _UM has been replaced by ΔT _UM ≈ \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\epsilon$$ \end{document} _UMΔT _m, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\epsilon$$ \end{document} _UM = exp(−(R _UM − R _m)αg/c_p ) ≈ 0.7 is the adiabatic temperature decrease from the average mantle temperature to the bottom of the upper mantle thermal boundary layer, ΔT _m = T _m − T_g , and the mantle cooling exponent is β = 1/3.

Radiogenic heat production in Earth is generated primarily by the decay of ²³⁸U, ²³⁵U, ²³²Th, and ⁴⁰K, which is approximated in the mantle by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}Q_{\rm rad} ( t ) = Q_{\rm{rad , 0}} \ { \rm exp} (- t / \tau_{\rm rad} ) \tag{26}\end{align*} \end{document}

where Q _rad,0 is the initial radiogenic heat production rate at t = 0 and τ _rad is the radioactive decay timescale that approximates the decay of the four major isotopes. The precise bulk silicate Earth radiogenic heat production rate is somewhat uncertain, so we use a nominal value of Q _rad(t =4.5 Gyr) = 13 TW (Jaupart et al., 2007).

Similar to the mantle convective heat flow, the CMB heat flow is \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}Q_ { \rm cmb } = A_ { \rm c } k_ { \rm LM } \frac { \Delta T_ { \rm LM } } { \delta_ { \rm LM } } \tag { 27 } \end{align*} \end{document}

where A _c is core surface area and k _LM is lower mantle thermal conductivity. The lower mantle and CMB temperatures, T _LM and T _cmb, are extrapolations along the mantle and core adiabats: T _LM = \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\epsilon$$ \end{document} _LM T _m and T _cmb = \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\epsilon$$ \end{document} _c T _c, where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\epsilon$$ \end{document} _LM = exp(−(R _LM − R _m)αg/c_p ) ≈ 1.3 and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\epsilon$$ \end{document} _c ≈ 0.8. The lower mantle thermal boundary layer thickness is also derived by assuming the boundary layer Rayleigh number is critical and that ν _LM = 2ν _UM, which was found by Driscoll and Bercovici (2014) to produce a nominal Earth model.

Core secular cooling is \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}Q_{\rm core} = - M_{\rm c} c_{\rm c} \dot {T}_{\rm c} \tag{28}\end{align*} \end{document}

where M _c is core mass, c _c is core specific heat, and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\dot {T}_{ \rm c}$$ \end{document} is the rate of change of the average core temperature T _c.

Radiogenic heat in the core is produced primarily by the decay of ⁴⁰K (Gessmann and Wood, 2002; Murthy et al., 2003; Corgne et al., 2007). Its time dependence is treated the same as mantle radiogenic heat in (26) but with a radioactive decay timescale of τ _rad,c = 1.2 Gyr. We assume an abundance of ⁴⁰K in the core that corresponds to 2 TW of heat production after 4.5 Gyr.

B.3. Melting

The mantle solidus is approximated by a third-order polynomial (Elkins-Tanton, 2008), \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}T_{\rm sol} ( r ) = A_{\rm sol} r^3 + B_{\rm sol} r^2 + C_{\rm sol} r + D_{\rm sol} \tag{29}\end{align*} \end{document}

where the coefficients are constants (see Table 1). This solidus is calibrated to fit the following constraints: solidus temperature of 1450 K at the surface, solidus temperature of 4150 K at the CMB (Andrault et al., 2011), and present-day upwelling melt fraction of f _melt = 8%. The liquidus is assumed to be hotter by a constant offset ΔT _liq = 500 K, so T _liq(r) = T _sol(r) + ΔT _liq.

Mantle melt heat loss (or advective heat flow) is modeled as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}Q_{\rm melt} = \epsilon_{\rm erupt} \dot {M}_{\rm melt} ( L_{\rm melt} + c_{\rm m} \Delta T_{\rm melt} ) \tag{30}\end{align*} \end{document}

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\epsilon_{ \rm erupt} = 0.2$$ \end{document} is the efficiency of magma eruption to the surface (assumed to be constant and equal to present-day value), \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\dot {M}_{ \rm melt}$$ \end{document} is melt mass flux (see below), L _melt is latent heat of the melt, c _m is specific heat of the melt, and ΔT _melt is the excess temperature of the melt at the surface (see below). This formulation of heat loss is similar to the “heat pipe” mechanism invoked for Io (O'Reilly and Davies, 1981; Moore, 2003), where melt is a significant source of heat loss. We note that this mechanism is more important for stagnant-lid planets where the normal conductive heat flow is lower (Driscoll and Bercovici, 2014).

The melt mass flux \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\dot {M}_{ \rm melt}$$ \end{document} is the product of the upwelling solid mass flux times the melt mass fraction f _melt, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\dot {M}_{\rm melt} = \dot {V}_{\rm up} \rho_{\rm solid} f_{\rm melt} ( z_{\rm UM} ) \tag{31}\end{align*} \end{document}

where solid density is ρ _solid, volumetric upwelling rate is \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\dot {V}$$ \end{document} _up = 1.16κA _p/δ _UM, z _UM = R − δ _UM, and melt fraction is \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}f_ { \rm melt } ( z ) = \frac { T_ { \rm m } ( z ) - T_ { \rm sol } } { T_ { \rm liq } - T_ { \rm sol } } \tag { 32 } \end{align*} \end{document}

This model predicts a ridge melt production of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\dot {M}$$ \end{document} _melt =2.4 × 10⁶ kg s⁻¹ for δ _UM = 80 km and f _melt = 0.1, similar to present-day global melt production estimates (Cogné and Humler, 2004).

We define the magma ocean as the region of the mantle with temperature exceeding the liquidus. Given the geotherm in (18, 20) and the liquidus T _liq(r) similar to (29), the mantle will mainly freeze from the bottom of the convecting mantle up because the liquidus gradient is steeper than the adiabat (e.g., Elkins-Tanton, 2012). However, if the core is hot enough, a second melt region exists in the lower mantle boundary layer, where the temperature gradient exceeds the liquidus and the mantle freezes toward the CMB. As can be seen in Fig. B1, a basal magma ocean exists for about 4 Gyr before solidifying.

FIG. B1.

Geotherm temperature profile T(r) for the orbit-evolving model with initial orbit of a = 0.05 AU and e = 0.5. Temperature profiles are shown at 1 Myr, 1 Gyr, 4.5 Gyr, and 10 Gyr. The core liquidus and mantle liquidus and solidus are shown in red.

Latent heat released from the solidification of the mantle is \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}Q_{\rm{L , man}} = \dot {M}_{\rm sol} L_{\rm melt} \tag{33}\end{align*} \end{document}

where L _melt is the latent heat released per kilogram and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\dot {M}_{ \rm sol}$$ \end{document} is the solid mantle growth rate. The growth rate is calculated assuming a uniform mantle density ρ _m so that \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\dot {M}_{ \rm sol} = \rho_{ \rm m} \dot {V}_{ \rm sol}$$ \end{document} , where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\dot {V}_{ \rm sol} = - \dot {V}_{ \rm liq}$$ \end{document} . The rate of change of the liquid volume of the mantle is \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\dot { V } _ { \rm liq } = \frac { dV_ { \rm liq } } { dT_ { \rm m } } \dot { T } _ { \rm m } \tag { 34 } \end{align*} \end{document}

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\dot {T}_{ \rm m}$$ \end{document} is the mantle secular cooling rate and dV _liq/dT _m is linearly approximated by 8 × 10¹⁷ m³ K⁻¹, which is the change in liquid volume from a 90% liquid to a completely solid mantle. This approximation implies that the latent heat released due to mantle solidification is linearly proportional to the mantle secular cooling rate, and the ratio of the latent heat flow to the mantle secular cooling heat flow is Q _L,man /Q _sec,m ≈ 0.24. For example, a mantle solidification time of 100 Myr corresponds to an average latent heat release of Q _L,man ≈ 400 TW over that time.

B.4. Core dynamo

Given the thermal cooling rate of the core, the magnetic dipole moment \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${\cal M}$$\end{document} \end{document} is estimated from the empirical scaling law, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}{ {\cal M}} = 4 \pi R_{\rm c}^3 \gamma_{\rm d} \sqrt{ \rho / 2 \mu_0} ( F_{\rm c} D_{\rm c} ) ^{1 / 3} \tag{35}\end{align*} \end{document}

where γ _d = 0.2 is the saturation constant for fast-rotating dipolar dynamos, μ ₀ = 4π × 10⁻⁷ H m⁻¹ is magnetic permeability, D _c = R _c − R _ic is the dynamo region shell thickness, R _c and R _ic are outer and inner core radii, respectively, and F _c is the core buoyancy flux (Olson and Christensen, 2006). We assume that the field is dipolar, ignoring the complicating influences of shell thickness and heterogeneous boundary conditions (e.g., Heimpel et al., 2005; Aubert et al., 2009; Driscoll and Olson, 2009; Olson et al., 2014). In this formulation, a positive buoyancy flux implies dynamo action, which is a reasonable approximation when the net buoyancy flux is large but may overestimate the field strength at low flux. The total core buoyancy flux F _c is the sum of thermal and compositional buoyancy fluxes, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}F_{\rm c} = F_{\rm th} + F_{\rm X} \tag{36}\end{align*} \end{document}

where the thermal and compositional buoyancy fluxes are \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}F_ {\rm th} = \frac { \alpha_{\rm c} q_ {\rm c}} {\rho_ {\rm c } c_ {\rm c } } q_{\rm { c , conv } }\tag { 37 } \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}F_ { \rm X } = g_ { \rm ic } \frac { \Delta \rho_X } { \rho_ { \rm c } } \left( \frac { R_ { \rm ic } } { R_ { \rm c } } \right) ^2 \ \dot { R } _ { \rm ic } \tag { 38 } \end{align*} \end{document}

where the subscript c refers to bulk core properties, core convective heat flux is q _c,conv = q _cmb − q _c,ad, gravity at the ICB is approximated by g _ic = g _c R _ic/R _c, and the outer core compositional density difference is Δρ_χ = ρ _c − ρ_χ with ρ_χ the light element density. For simplicity, the expression for light element buoyancy (38) ignores buoyancy due to latent heat release at the ICB because it is a factor of 3.8 less than buoyancy of the light elements.

The isentropic core heat flux at the CMB, proportional to the gradient of (22), is \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}q_{\rm{c , ad}} = k_{\rm c} T_{\rm cmb}R_{\rm c} / D^2_N \tag{39}\end{align*} \end{document}

where core thermal conductivity is approximated by the Wiedemann-Franz law, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}k_{\rm c} = \sigma_{\rm c} L_{\rm c} T_{\rm cmb} \tag{40}\end{align*} \end{document}

and electrical conductivity is σ _c and L _c is the Lorentz number. For typical values of high pressure-temperature iron, σ _c = 10 × 10⁵ Ω⁻¹ m⁻¹ (Pozzo et al., 2012; Gomi et al., 2013), L _c = 2.5 × 10⁻⁸ W Ω K⁻¹, and T _cmb = 4000 K, the core thermal conductivity is k _c = 100 W m⁻¹ K⁻¹.

Footnotes

Acknowledgments

The authors thank W. Henning and T. Hurford for helpful discussions and two anonymous reviewers for valuable feedback. This work was performed as part of the NASA Astrobiology Institute's Virtual Planetary Laboratory, supported by the National Aeronautics and Space Administration through the NASA Astrobiology Institute under solicitation NNH12ZDA002C and Cooperative Agreement Number NNA13AA93A.

Author Disclosure Statement

No competing financial interests exist.

Abbreviations Used

References

Andrault

, Bolfan-Casanova

, Nigro

G.L.

, Bouhifd

M.A.

, Garbarino

, and Misra

(2011) Solidus and liquidus profiles of chondritic mantle: implication for melting of the Earth across its history. Earth Planet Sci Lett, 304:251–259.

Aubert

, Labrosse

, and Poitou

(2009) Modeling the paleo-evolution of the geodynamo. Geophys J Int, 179:1414–1428.

Barnes

and Greenberg

(2006) Stability limits in extrasolar planetary systems. Astrophys J, 647, doi:10.1086/507521.

Barnes

and Heller

(2013) Habitable planets around white and brown dwarfs: the perils of a cooling primary. Astrobiology, 13:279–291.

Barnes

and Quinn

(2004) The (in)stability of planetary systems. Astrophys J, 611:494–516.

Barnes

, Raymond

S.N.

, Greenberg

, Jackson

, and Kaib

N.A.

(2010) CoRoT-7b: super-Earth or super-Io?. Astrophys J, 709, doi:10.1088/2041-8205/709/2/L95.

Barnes

, Mullins

, Goldblatt

, Meadows

V.S.

, Kasting

J.F.

, and Heller

(2013) Tidal Venuses: triggering a climate catastrophe via tidal heating. Astrobiology, 13:225–250.

Běhounková

, Tobie

, Choblet

, and Čadek

(2010) Coupling mantle convection and tidal dissipation: applications to Enceladus and Earth-like planets. J Geophys Res Planets, 115, doi:10.1029/2009JE003564.

Běhounková

, Tobie

, Choblet

, and Čadek

(2011) Tidally induced thermal runaways on extrasolar Earths: impact on habitability. Astrophys J, 728, doi:10.1088/0004-637X/728/2/89.

10.

Cogné

J.-P.

and Humler

(2004) Temporal variation of oceanic spreading and crustal production rates during the last 180 My. Earth Planet Sci Lett, 227:427–439.

11.

Corgne

, Keshav

, Fei

, and McDonough

W.F.

(2007) How much potassium is in the Earth's core? New insights from partitioning experiments. Earth Planet Sci Lett, 256:567–576.

12.

Costa

, Caricchi

, and Bagdassarov

(2009) A model for the rheology of particle-bearing suspensions and partially molten rocks. Geochemistry, Geophysics, Geosystems, 10, doi:10.1029/2008GC002138.

13.

Dartnell

L.R.

(2011. Ionizing radiation and life. Astrobiology, 11:551–582.

14.

Dole

S.H.

(1964) Habitable Planets for Man, Blaisdell, New York.

15.

Driscoll

and Bercovici

(2013) Divergent evolution of Earth and Venus: influence of degassing, tectonics, and magnetic fields. Icarus, 226:1447–1464.

16.

Driscoll

and Bercovici

(2014) On the thermal and magnetic histories of Earth and Venus: influences of melting, radioactivity, and conductivity. Physics of the Earth and Planetary Interiors, 236:36–51.

17.

Driscoll

and Olson

(2009) Effects of buoyancy and rotation on the polarity reversal frequency of gravitationally driven numerical dynamos. Geophys J Int, 178:1337–1350.

18.

Driscoll

and Olson

(2011) Optimal dynamos in the cores of terrestrial exoplanets: magnetic field generation and detectability. Icarus, 213:12–23.

19.

Dwyer

, Stevenson

, and Nimmo

(2011) A long-lived lunar dynamo driven by continuous mechanical stirring. Nature, 479:212–214.

20.

Elkins-Tanton

L.T.

(2008) Linked magma ocean solidification and atmospheric growth for Earth and Mars. Earth Planet Sci Lett, 271:181–191.

21.

Elkins-Tanton

L.T.

(2012) Magma oceans in the inner solar system. Annu Rev Earth Planet Sci, 40:113–139.

22.

Ferraz-Mello

, Rodríguez

, and Hussmann

(2008) Tidal friction in close-in satellites and exoplanets: the Darwin theory re-visited. Celestial Mechanics and Dynamical Astronomy, 101:171–201.

23.

Foley

B.J.

, Bercovici

, and Landuyt

(2012) The conditions for plate tectonics on super-Earths: inferences from convection models with damage. Earth Planet Sci Lett, 331–332:281–290.

24.

Gaidos

, Conrad

, Manga

, and Hernlund

(2010) Thermodynamic limits on magnetodynamos in rocky exoplanets. Astrophys J, 718:596–609.

25.

Gessmann

and Wood

(2002) Potassium in the Earth's core?. Earth Planet Sci Lett, 200:63–78.

26.

Goldblatt

and Watson

A.J.

(2012) The runaway greenhouse: implications for future climate change, geoengineering and planetary atmospheres. Philos Transact A Math Phys Eng Sci, 370:4197–4216.

27.

Goldreich

and Soter

(1966) Q in the Solar System. Icarus, 5:375–389.

28.

Gomi

, Ohta

, Hirose

, Labrosse

, Caracas

, Verstraete

M.J.

, and Hernlund

J.W.

(2013) The high conductivity of iron and thermal evolution of the Earth's core. Physics of the Earth and Planetary Interiors, 224:88–103.

29.

Griessmeier

J.-M.

, Stadelmann

, Motschmann

, Belisheva

, Lammer

, and Biernat

(2005) Cosmic ray impact on extrasolar Earth-like planets in close-in habitable zones. Astrobiology, 5:587–603.

30.

Heimpel

, Aurnou

, Al-Shamali

, and Gomez Perez

(2005) A numerical study of dynamo action as a function of spherical shell geometry. Earth Planet Sci Lett, 236:542–557.

31.

Henning

W.G.

and Hurford

(2014) Tidal heating in multilayered terrestrial exoplanets. Astrophys J, 789, doi:10.1088/0004-637X/789/1/30.

32.

Henning

W.G.

, O'Connell

, and Sasselov

(2009) Tidally heated terrestrial exoplanets: viscoelastic response models. Astrophys J, 707:1000–1015.

33.

Hirose

, Labrosse

, and Hernlund

(2013) Composition and state of the core. Annu Rev Earth Planet Sci, 41:657–691.

34.

Howard

L.N.

(1966) Convection at high Rayleigh number. In Applied Mechanics, Springer, Berlin, pp 1109–1115.

35.

Jackson

, Barnes

, and Greenberg

(2009) Observational evidence for tidal destrucion of exoplanets. Astrophys J, 698:1357–1366.

36.

Jackson

, Faul

U.H.

, Gerald

, John

, and Tan

B.H

. (2004) Shear wave attenuation and dispersion in melt-bearing olivine polycrystals: 1. Specimen fabrication and mechanical testing. J Geophys Res Solid Earth, 109, doi:10.1029/2003JB002406.

37.

Jaupart

, Labrosse

, and Mareschal

J.C.

(2007) Temperatures, heat and energy in the mantle of the Earth. In Treatise on Geophysics, Vol. 7: Mantle Dynamics, edited by Bercovici

, Elsevier, Amsterdam, pp 253–305.

38.

Kaula

W.M.

(1964) Tidal dissipation by solid friction and the resulting orbital evolution. Rev Geophys, 2:661–685.

39.

Kopparapu

R.K.

, Ramirez

, Kasting

J.F.

, Eymet

, Robinson

T.D.

, Mahadevan

, Terrien

R.C.

, Domagal-Goldman

, Meadows

, and Deshpande

(2013) Habitable zones around main-sequence stars: new estimates. Astrophys J, 765, doi:10.1088/0004-637X/765/2/131.

40.

Kopparapu

R.K.

, Ramirez

R.M.

, SchottelKotte

, Kasting

J.F.

, Domagal-Goldman

, and Eymet

(2014) Habitable zones around main-sequence stars: dependence on planetary mass. Astrophys J, 787, doi:10.1088/2041-8205/787/2/L29.

41.

Korenaga

(2006) Archean geodynamics and the thermal evolution of Earth. In Archean Geodynamics and Environments, Geophysical Monograph No. 164, edited by Benn

, Mareschal

J.-C.

, and Condie

K.C.

, American Geophysical Union, Washington, DC, pp 7–32.

42.

Labrosse

, Poirier

, and Le Mouël

(2001) The age of the inner core. Earth Planet Sci Lett, 190:111–123.

43.

Le Bars

, Cébron

, and Le Gal

(2015) Flows driven by libration, precession, and tides. Annu Rev Fluid Mech, 47:163–193.

44.

Luger

, Barnes

, Lopez

, Fortney

, Jackson

, and Meadows

(2015) Habitable evaporated cores: transforming mini-Neptunes into super-Earths in the habitable zones of M dwarfs. Astrobiology, 15:57–88.

45.

Matsuyama

(2014) Tidal dissipation in the oceans of icy satellites. Icarus, 242:11–18.

46.

Mayor

, Lovis

, and Santos

N.C.

(2014) Doppler spectroscopy as a path to the detection of Earth-like planets. Nature, 513:328–335.

47.

McWilliams

J.C.

(2012) The elemental shear dynamo. J Fluid Mech, 699:414–452.

48.

Misra

, Krissansen-Totton

, Koehler

M.C.

, and Sholes

(2015) Transient sulfate aerosols as a signature of exoplanet volcanism. arXiv:1504.04629.

49.

Moore

W.B.

(2003) Tidal heating and convection in Io. J Geophys Res Planets, 108, doi:10.1029/2002JE001943.

50.

Moore

W.B.

and Webb

A.A.G.

(2013) Heat-pipe Earth. Nature, 501:501–505.

51.

Moore

W.B.

, Schubert

, Anderson

J.D.

, and Spencer

J.R.

(2007) The interior of Io. In Io After Galileo, edited by Lopes

R.M.C.

and Spencer

J.R.

, Springer, Berlin, pp 89–108.

52.

Murthy

V.R.

, van Westrenen

, and Fei

(2003) Experimental evidence that potassium is a substantial radioactive heat source in planetary cores. Nature, 423:163–165.

53.

Olson

and Christensen

(2006) Dipole moment scaling for convection-driven planetary dynamos. Earth Planet Sci Lett, 250:561–571.

54.

Olson

, Deguen

, Hinnov

, and Zhong

(2013) Controls on geomagnetic reversals and core evolution by mantle convection in the Phanerozoic. Physics of the Earth and Planetary Interiors, 214:87–103.

55.

Olson

, Hinnov

L.A.

, and Driscoll

P.E.

(2014) Nonrandom geomagnetic reversal times and geodynamo evolution. Earth Planet Sci Lett, 388:9–17.

56.

O'Reilly

T.C.

and Davies

G.F.

(1981) Magma transport of heat on Io: a mechanism allowing a thick lithosphere. Geophys Res Lett, 8:313–316.

57.

Peale

and Cassen

(1978) Contribution of tidal dissipation to lunar thermal history. Icarus, 36:245–269.

58.

Pozzo

, Davies

, Gubbins

, and Alfè

(2012) Thermal and electrical conductivity of iron at Earth's core conditions. Nature, 485:355–358.

59.

Rasio

, Tout

, Lubow

, and Livio

(1996) Tidal decay of close planetary orbits. Astrophys J, 470:1187–1191.

60.

Ray

R.D.

and Egbert

G.D.

(2012) Fortnightly Earth rotation, ocean tides and mantle anelasticity. Geophys J Int, 189:400–413.

61.

Reichow

M.K.

, Saunders

A.D.

, White

R.V.

, Pringle

M.S.

, Al'Mukhamedov

A.I.

, Medvedev

A.I.

, and Kirda

N.P

. (2002) ⁴⁰Ar/³⁹Ar dates from the West Siberian Basin: Siberian flood basalt province doubled. Science, 296:1846–1849.

62.

Seager

(2013) Exoplanet habitability. Science, 340:577–581.

63.

Segatz

, Spohn

, Ross

, and Schubert

(1988) Tidal dissipation, surface heat flow, and figure of viscoelastic models of Io. Icarus, 75:187–206.

64.

Solomatov

(1995) Scaling of temperature-and stress-dependent viscosity convection. Physics of Fluids, 7, doi:10.1063/1.868624.

65.

Solomatov

and Moresi

(2000) Scaling of time-dependent stagnant lid convection: application to small-scale convection on Earth and other terrestrial planets. J Geophys Res, 105:21795–21817.

66.

Sotin

and Labrosse

(1999) Three-dimensional thermal convection in an iso-viscous, infinite Prandtl number fluid heated from within and from below: applications to the transfer of heat through planetary mantles. Physics of the Earth and Planetary Interiors, 112:171–190.

67.

Sotin

, Tobie

, Wahr

, and McKinnon

W.B.

(2009) Tides and tidal heating on Europa. In Europa, edited by Pappalardo

R.T.

, McKinnon

W.B.

, and Khurana

, University of Arizona Press, Tucson, pp 85–117.

68.

Tachinami

, Senshu

, and Ida

(2011) Thermal evolution and lifetime of intrinsic magnetic fields of super-Earths in habitable zones. Astrophys J, 726, doi:10.1088/0004-637X/726/2/70.

69.

Tyler

(2014) Comparative estimates of the heat generated by ocean tides on icy satellites in the outer Solar System. Icarus, 243:358–385.

70.

Van Laerhoven

, Barnes

, and Greenberg

(2014) Tides, planetary companions, and habitability: habitability in the habitable zone of low-mass stars. Mon Not R Astron Soc, 441:1888–1898.

71.

Van Summeren

, Gaidos

, and Conrad

C.P.

(2013) Magnetodynamo lifetimes for rocky, Earth-mass exoplanets with contrasting mantle convection regimes. J Geophys Res Planets, 118:938–951.

72.

Zimmerman

D.S.

, Triana

S.A.

, Nataf

H.-C.

, and Lathrop

D.P

. (2014) A turbulent, high magnetic Reynolds number experimental model of Earth's core. J Geophys Res Solid Earth, 119:4538–4557.

73.

Zuluaga

J.I.

, Bustamante

, Cuartas

P.A.

, and Hoyos

J.H.

(2013) The influence of thermal evolution in the magnetic protection of terrestrial planets. Astrophys J, 770, doi:10.1088/0004-637X/770/1/23.