THERMODYNAMIC LIMITS ON MAGNETODYNAMOS IN ROCKY EXOPLANETS

A planet is modeled as a homogeneous, fully convecting Mg/Fe–silicate mantle surrounding a liquid/solid Fe core. The pressure–density profiles are calculated self-consistently with the (time-varying) size of the inner core. Third-order Birch–Murnagham (BM) equations of state (EOS) are used for each component. The pressure is

$$\begin{equation} P = \frac{3}{2}K_1(x^7-x^5)\left[1 + \frac{3}{4}(4-K_1^{\prime })(1-x^2)\right], \end{equation} \tag{ 6 }$$

where x = (ρ/ρ₁)^1/3 is the isotropic strain. Values of the ambient density ρ₁, compressibility K₁, and its pressure-derivative K'₁ for the relevant planetary materials are given in Table 1. Electron degeneracy (Thomas–Fermi–Dirac, or TFD) pressure will be important in the metal cores of massive Earth-like planets. We calculated the TFD contribution to the pressure using the formulation of Zapolsky & Salpeter (1969) and found that the EOS for liquid and -phase solid iron (Fe(l) and Fe()) intersect at about 3.3 TPa. Pressures beyond 3 TPa are not relevant to the range of masses considered here. The density in the liquid part of the core is adjusted by a fixed fraction δρ/ρ to account for the presence of light elements (Li & Fei 2003).

Table 1. EOS Parameters

Material	ρ (kg m−3)	K1 (GPa)	K'1	Reference
Ice/water	1000	2.2	4.0	Sotin et al. (2007)
HP ice (VII)	1460	23.9	4.2	Sotin et al. (2007)
Mg–Olivine	3222	128	4.3	Sotin et al. (2007)
Perovskite	4260	266	3.9	Ahrens (2000)
Fe(l)	7019	154a	4.66	Anderson & Ahrens (1994)
Fe()	8315	165	4.97	Vocadlo et al. (2003)

Note. ^aAdjusted 40% upward from the Anderson & Ahrens (1994) value.

Download table as:  ASCIITypeset image

Our interior model reproduces the central pressure (364 GPa), CMB pressure (136 GPa), and core size R_c (3480 km) of the Earth (Dziewonski & Anderson 1981), but only if the zero-pressure compressibility for Fe(l) is 40% larger than the 106 GPa value reported by Anderson & Ahrens (1994) based on shock experiments. With this adjustment, the predicted radius of the Earth is 25 km (0.4%) less than the actual value, but this is not surprising because we do not include low-density crustal rocks or an ocean. Our dynamo predictions are much more sensitive to the properties of the core than those of the entire planet. Predicted planet radius, core radius, pressure at the CMB, and central pressure are plotted versus total mass in Figure 2. Pressures increase as the core solidifies and contracts.

Zoom In Zoom Out Reset image size

The adiabatic temperature profile in the outer convecting core is approximately

$$\begin{equation} T(r) \approx T_c \exp \left[(R_c^2 - r^2)/D_0^2\right]. \end{equation} \tag{ 7 }$$

The thermal length scale evaluated at the planet center is

$$\begin{equation} D_0 = \sqrt{\frac{3c_p}{2\pi \alpha _0 \rho _0 G}}, \end{equation} \tag{ 8 }$$

(Labrosse 2003), where c_p and α are the heat capacity and thermal expansivity of Fe, respectively, and G is the gravitational constant (see Table 2 for parameter values). The temperature profile deviates from Equation (7) to the extent that D, and specifically α, vary with depth in the outer core. α must therefore be evaluated appropriately when calculating terms that are sensitive to D, i.e., the heat and entropy transported by thermal conduction. Hereafter, D (without a subscript) is D₀.

The density in the core follows

$$\begin{equation} \rho (r) \approx \rho _c \exp \left[\big(R_c^2 - r^2\big)/L^2\right], \end{equation} \tag{ 9 }$$

where the density scale length is

$$\begin{equation} L = \sqrt{\frac{9K_1}{2\pi G \rho _1 \rho _0} \left(\ln \frac{\rho _0}{\rho _1} + 1 \right)}. \end{equation} \tag{ 10 }$$

For the Earth, D ∼ 6400 km and L ∼ 7400 km and these values are only weakly dependent on planet size. We expect the core radius in planets of Earth-like composition to scale as M^0.25_p (Sotin et al. 2007; Seager et al. 2007). The nth order term in either r/D or r/L will be of order M^n/4_p/2ⁿ. For planets of a few Earth masses, we must retain terms in expansions with r/D or r/L through fifth order.

If a solid inner core is present, the temperature at the ICB is the intercept between the adiabat and the iron melting point τ. Adopting a Lindemann law for τ,

$$\begin{equation} \frac{\partial \log \tau }{\partial \log \rho } = 2\left(\gamma - \frac{1}{3}\right), \end{equation} \tag{ 11 }$$

where γ is the Grüneisen parameter. The application of the Lindemann law to the deep interiors of planets is sometimes viewed as speculative (e.g., Wolf & Jeanloz 1984). The assumptions upon which it is derived are not strictly valid for polyatomic systems, and the law is known to sometimes fail for minerals with complex interatomic forces and structures. However, at high pressures and for metals it becomes a better approximation, and hence is widely used to interpret shock data (e.g., Anderson & Ahrens 1996). Importantly, it appears to provide a good description of experimental data and first principle calculations for Fe (e.g., Hemley & Mao 2001; Huser et al. 2005) and other metals (e.g., Dai et al. 2002).

The radial dependence of the solidus becomes

$$\begin{equation} \tau (r) = \tau _0 \exp \left[-2\left(1-\frac{1}{3\gamma }\right)\left(\frac{r}{D}\right)^2\right], \end{equation} \tag{ 12 }$$

where we have used L²/D² = γ (Labrosse 2003). Thus, the temperature at the top of the core is uniquely set by the ratio of the inner to outer core radius ℜ = R_i/R_c:

$$\begin{equation} T_c = \tau _0 \exp \left[-\left[1 + \Re ^2 \left(1 - \frac{2}{3\gamma }\right)\right]\frac{R_c^2}{D^2}\right]. \end{equation} \tag{ 13 }$$

T_c must decrease as the inner core grows. If the entire core is liquid, then T_c is the independent variable. For greater accuracy, we use our interior model to directly calculate the density at the ICB, the melting temperature using Equation (11), and T_c using Equation (7).

The temperature profile of the inner core will depend on whether it, too, is convecting (Buffett 2009). Such convection may be sufficiently torpid so as not to contribute to the dynamo, but as long as overturn is rapid compared to the timescale of inner core growth (10⁹–10¹⁰ yr), the temperature profile will be an adiabat. This assumption maximizes the transport of sensible heat into the outer core, and hence minimizes the rate of inner core growth.

For a fully convecting core, the heat conducted along the adiabat at the CMB is (neglecting compressibility effects)

$$\begin{equation} Q_K \approx \frac{4 \pi R_c^2 k \alpha _c g_c T_c}{c_p}, \end{equation} \tag{ 14 }$$

where k is the thermal conductivity and g is the local gravitational acceleration. If Q_c, the heat flow from the core to the mantle across the CMB, is less than this, the temperature profile will be subadiabatic to the degree required to equalize the heat flows. Whether this layer is completely stratified will depend on non-thermal (i.e., light element) buoyancy forces. In the absence of a growing inner core, however, such a layer will be thermally stratified. Geomagnetic variation limits the thickness of any upper stratified layer in the Earth's core to <100 km (Gubbins 2007), if one exists at all (Stanley & Mohammadi 2008).

If the stratified layer is thin and makes a negligible contribution to the sensible heat flow, the heat flow through it is constant and its temperature profile is simple. The location ℜ_* and temperature T_* of the boundary between the convective and conducting regions is found by simultaneously solving for the heat flow and temperature:

$$\begin{equation} \Re _*^3 = \left(\frac{D_c}{R_c}\right)^3\frac{Q_c(T_c)}{8 \pi k D_c T_*} \end{equation} \tag{ 15 }$$

and

$$\begin{equation} T_c = T_* - \frac{Q_c (T_c)}{4\pi k R_c}\left(\frac{1}{\Re _*}-1\right), \end{equation} \tag{ 16 }$$

and using Equation (7), i.e.,

$$\begin{equation} T_* = T_i \exp \left[\left(\Re ^2-\Re _*^2\right)\left(\frac{R_c}{D_c}\right)^2\right], \end{equation} \tag{ 17 }$$

where the dependence of Q_c on T_c must be accounted for (Section 3.4), and we evaluate the thermal scale length at the CMB. In the event that stratification does occur, we compute the entropy terms for the convecting part of the core only. This is done by re-scaling the radius and mass of the core in the entropy equations to the size and mass of the convecting zone. T_c becomes the temperature T_* at the top of the convecting zone, and we assume that since the adiabat is self-similar, such that dT_*/T_* ≈ dT_c/T_c (effectively ignoring small differences due to changes in the thickness of the convective zone itself).

The condition for maintenance of convective motions for a dynamo can be expressed as a balance between sources and sinks of entropy in the core, e.g., Labrosse (2007) and Nimmo (2009):

$$\begin{equation} \Phi + E_{K} = E_R + E_S + E_G + E_L, \end{equation} \tag{ 18 }$$

where E_K is the entropy sink due to conduction along the adiabat, E_R the production of entropy by internal heating, E_S that due to core cooling (loss of sensible heat), E_G the entropy production from buoyancy generation, and E_L the entropy production due to latent heat generated by crystallization of the core. The decay of radioisotopes is a possible internal heat source. ⁶⁰Fe will decay completely (t_1/2 = 2.6 Myr; Rugel 2009) before planets finish accreting, but ⁴⁰K (t_1/2 = 1.25 Gyr) may be a heat source in Earth's core (Nimmo et al. 2004). Values of the partition coefficient of K between silicate and Fe/Ni alloy melts have been discrepant between different high-pressure experiments, possibly for technical reasons (Li & Fei 2003), but new results limit the abundance of K to a few tens of parts per million, with no evidence for an increase in solubility with pressure (Corgne et al. 2007). In the absence of radiogenic heating,

$$\begin{equation} \Phi = E_L + E_S + E_G - E_K. \end{equation} \tag{ 19 }$$

We calculate Φ, divide by the outer core mass to obtain ϕ, and use Equations (2), (4), and (5) to estimate the magnetic field intensity and magnetic Reynolds number.

Expressions for entropy production (dissipation) have been derived and discussed extensively in the literature, e.g., Lister (2003). Following Labrosse (2007) and Nimmo (2009), we derive expansions in terms of the size of the inner core. The term from the loss of sensible heat comes form the thermodynamic identity dS = c_pdT/T under constant pressure and no heats of reaction. Thus,

$$\begin{equation} E_S = -\int \left(\frac{1}{T_c} - \frac{1}{T}\right) \rho c_p \frac{\partial T}{\partial t} dV, \end{equation} \tag{ 20 }$$

where the integration is over the volume of the entire core. In the absence of an inner core, and to fifth order in R_c/D and R_c/L,

$$\begin{equation} E_S \approx - \frac{2}{5}c_p M_c \left(\frac{R_c}{D}\right)^2 \left[1 + \frac{2}{7}\left(\frac{R_c}{D}\right)^2 + \frac{6}{35}\left(\frac{R_c}{L}\right)^2\right]\frac{dT_c}{dt}, \end{equation} \tag{ 21 }$$

where M_c is the total mass of the core. We assume that subsolidus convection continues the adiabatic temperature profile into the inner core (Buffett 2009) and ignore the effect of small changes in the temperature and pressure length scales due to differences in compressibility and density between the liquid and solid () phases. This assumption maximizes the contribution of the inner core to the entropy and energy terms. In the case of a mostly solid core, a detailed treatment of the heat flow from the inner core would be warranted. However, we shall show that in planets much larger than the Earth, inner cores remain relatively small, or do not grow at all.

The entropic term for the latent heat due to core growth is

$$\begin{equation} E_L = 4\pi R_i^2 \rho _i \Delta S \left(\frac{T_i}{T_c} - 1\right)\frac{dR_i}{dt}, \end{equation} \tag{ 22 }$$

where ΔS is the specific entropy of fusion. Assuming that the core grows outward,

$$\begin{equation} \frac{dR_i}{dt} = -\frac{D^2}{2R_i(\Delta -1)}\frac{1}{T_c}\frac{dT_c}{dt}, \end{equation} \tag{ 23 }$$

where Δ is the ratio of the temperature–pressure slopes of the solidus to the adiabat, evaluated at R_i. Then

$$\begin{eqnarray} &&E_L = -\frac{3}{2}\frac{M_c \Delta S }{\Delta - 1}\frac{\rho _i}{\bar{\rho }} \Re (1 - \Re ^2) \left[1 + \left(\frac{R_c}{D}\right)^2 \frac{1-\Re ^2}{2}\right]\frac{1}{T_c}\frac{dT_c}{dt},\nonumber \\ \end{eqnarray} \tag{ 24 }$$

where the ratio of the density at the top of the inner core to the mean density is obtained directly from the interior model (Section 3.1).

The entropic term for the release of buoyant fluid (density deficit Δρ) during the crystallization of the core is (Labrosse 2007)

$$\begin{equation} E_G \approx -\frac{3 \pi G M_c \Delta \rho D^2}{\Delta - 1}\frac{\rho _i}{\bar{\rho }}\Re \frac{F}{F^{\prime}}T_c^2\frac{dT_c}{dt}, \end{equation} \tag{ 25 }$$

where Δρ is the density change in the fluid due to the presence of light elements, and F and F' are the dimensionless functions

$$\begin{eqnarray} &&\begin{array}{c} \hspace{-7pc}\displaystyle{F = \frac{1}{5} - \frac{1}{3}\Re ^2 + \frac{2}{15}\Re ^{5}}\\[5pt] \displaystyle{+\left(\frac{R_c}{L} \right)^2 \left[\frac{2(1-\Re ^7)}{21} - \frac{1-\Re ^4}{6}\right]} \\[5pt] \displaystyle{+\left(\frac{R_c}{L}\right)^4 \left[\frac{1-\Re ^9}{27} - \frac{1-\Re ^6}{18} \right],} \\[5pt] \displaystyle{F^{\prime } = 1 \,{-}\, \Re ^3 \,{-}\, \frac{3}{5} \left(\frac{R_c}{L}\right)^2(1-\Re ^5) - \frac{3}{14}\left(\frac{R_c}{L}\right)^4(1-\Re ^7),} \end{array} \end{eqnarray} \tag{ 26 }$$

cf. Nimmo (2009). To derive this, we ignore work done by expansion of the core (Lister 2003), and assume that the energy is dissipated at temperature T_c.

The entropy production by the conduction of heat along the adiabat is

$$\begin{equation} E_K = \int _{0}^{R_c} 4\pi r^2 k \left(\frac{1}{T}\frac{\partial T}{\partial r}\right)^2 dr \end{equation} \tag{ 27 }$$

(Lister 2003). This is evaluated to be

$$\begin{equation} E_K = \frac{12 k M_c}{5\bar{\rho }D_c^2}\left(\frac{R_c}{D_c}\right)^2(1-\Re ^5). \end{equation} \tag{ 28 }$$

To more accurately reflect the temperature profile, to which conduction is sensitive, we have evaluated the temperature scale length at the CMB, rather than the core center. The other terms in the entropy equation depend on D to a lesser degree and are evaluated over the entire core, not just the CMB.

We solve for the core cooling rate dT_c/dt by balancing the heat flow:

$$\begin{equation} Q_S + Q_L + Q_G = Q_c, \end{equation} \tag{ 29 }$$

where the terms correspond to sensible heat, latent heat, gravitational energy, and heat transported through the CMB, respectively. The left-hand terms are all proportional to dT_c/dt. The sensible heat term is

$$\begin{equation} Q_S = \int \rho c_p \frac{\partial T}{\partial t} dV. \end{equation} \tag{ 30 }$$

Assuming the adiabat continues into the inner core, this is

$$\begin{equation} Q_S \approx c_p M_c \left[\! 1 + \frac{2}{5}\left(\frac{R_c}{D}\right)^2 \! + \frac{4}{35}\left(\frac{R_c}{D}\right)^4 \! + \frac{12}{175}\left(\frac{R_c^2}{DL}\right)^2\right]\!\frac{dT_c}{dt}. \end{equation} \tag{ 31 }$$

The latent heat term is

$$\begin{equation} Q_L = \frac{3}{2}\frac{M_c \Delta S \Re }{\Delta - 1} \left(\frac{D}{R_c}\right)^2 \exp \left[(R_c/D)^2(1-\Re ^2)\right]\frac{dT_c}{dt}, \end{equation} \tag{ 32 }$$

while the gravitational term is

$$\begin{equation} Q_G \approx -\frac{3 \pi G M_c \Delta \rho D^2}{\Delta - 1}\frac{\rho _i}{\bar{\rho }}\frac{F}{T_c}\frac{dT_c}{dt}. \end{equation} \tag{ 33 }$$

We compare Q_c to the heat Q_K that can be carried by conduction along the core adiabat:

$$\begin{equation} Q_K = \frac{8\pi R_c^3 k T_c}{D_c^2}. \end{equation} \tag{ 34 }$$

If Q_c > Q_K, then the entire core has thermally driven convection. If Q_c < Q_K, then cooling is insufficient to maintain convection in the entire core and a stratified, conducting layer develops (Section 3.2). In our calculation of the rate of sensible heat loss, we have neglected the small corrections due to departures of the temperature profile from the purely adiabatic case.

We assume that the mantle boundary layer at the CMB develops independently of the flow in the rest of the mantle and its thickness is controlled by a local Rayleigh number Ra, as opposed to the whole-mantle Ra (Nimmo & Stevenson 2000). Heat is removed from the core according to

$$\begin{equation} Q_c = 4\pi R_c k_m {\rm Nu}_c \left(T_c - T_l\right), \end{equation} \tag{ 35 }$$

where T_l is the temperature of the lower mantle and k_m is the thermal conductivity of the mantle. The Nusselt number is

$$\begin{equation} {\rm Nu}_c \approx ({\rm Ra}_c/{\rm Ra}_*)^{\delta }, \end{equation} \tag{ 36 }$$

where δ ≈ 0.3 (Schubert et al. 2001). For a fluid heated from below,

$$\begin{equation} {\rm Ra}_{c} = \frac{\rho g \alpha \left(T_c - T_l\right) d^3}{\kappa \eta _c}, \end{equation} \tag{ 37 }$$

where the gravitational acceleration g and viscosity η are evaluated at the CMB, and d = R_p − R_c. In the case of a partially stratified layer, we substitute T_t for T_c, where T_t is the temperature at the top of the stratified, conducting zone.

The heat transport by subsolidus convection in the mantle is

$$\begin{equation} Q_{m} = \frac{4 \pi R_p^2 {\rm Nu}_m k (T_m - T_p)}{R_p - R_c}, \end{equation} \tag{ 38 }$$

where T_p is the surface temperature. For a uniform internal heat source,

$$\begin{equation} {\rm Ra_m} = \frac{\rho g H \alpha d^5}{k \kappa \eta _m}, \end{equation} \tag{ 39 }$$

where H is the total specific heat generation (per unit mass). η_m is evaluated at the mantle reference temperature, taken to be halfway along the adiabat. We make the approximation that heat from the CMB is mixed rapidly and uniformly into the mantle (e.g., by plumes) and thus H = H_c + H_r, where H_c = Q_c/M_m, where M_m is the mass of the mantle, and the radiogenic heat production H_r is given by

$$\begin{equation} H_r = \Sigma _i C_i h_i e^{-t/\bar{t}_i}, \end{equation} \tag{ 40 }$$

where C_i is the initial concentration of the ith radioactive isotope, h_i is the specific heat production, and $\bar{t}_i$ is the mean life. The concentrations adopted for the long-lived radioisotopes ⁴⁰K, ²³²Th, ²³⁵U, and ²³⁸U are for the case of an "undepleted" terrestrial mantle (Ringwood 1991) given in Table 2 of Kite et al. (2009). The contribution due to mantle cooling is

$$\begin{equation} Q_T = -c_m M_m \frac{dT_m}{dt}, \end{equation} \tag{ 41 }$$

where c_m is the is the specific heat capacity of the mantle.

We adopt the viscosity law for a stress exponent n = 1 material, appropriate for the plate-tectonic regime (Nimmo & Stevenson 2000):

$$\begin{equation} \eta _m = \eta _* \exp \left(b \tau _m /T\right), \end{equation} \tag{ 42 }$$

where b ≈ 17 (Karato et al. 2001), τ_m is the peridotite solidus, and η_* is a reference viscosity. In Equation (37), we evaluate η at the "film" temperature, which is intermediate of the lower mantle and top of the core (Manga et al. 2001).

Convection in the silicate mantle may involve a mobile lid (i.e., plate tectonics) or a stagnant lid, and in our solar system we have Earth-mass examples of each (Earth and Venus). Whether one mode or another will be more likely on larger planets will depend on the relative magnitudes of lithosphere stresses (Valencia et al. 2007), the availability of water to weaken faults (Nimmo & McKenzie 1998), surface temperatures (Lenardic et al. 2008), and the formation of buoyant, unsubductable crust (Kite et al. 2009). In our model, the mode of mantle convection controls (1) the temperature boundary condition at the mantle side of the CMB and (2) the Nusselt–Rayleigh number relationship in the mantle.

In the case of a mobile lid, we assume that subducted slabs reach the lower mantle, have a temperature equal to the mean mantle temperature T_m, and directly cool the core. Thus, for Equation (35), we set T_l = T_m. We the use the same Nu–Ra relation as in Equation (36). In the stagnant lid case, the lower mantle temperature is set to

$$\begin{equation} T_l = \theta \exp \left[\int _{R_c}^{R_p} \frac{\alpha g}{c_p} dr\right], \end{equation} \tag{ 43 }$$

where θ is the potential temperature of the mantle, and the relationship between mean and potential temperatures is taken to be

$$\begin{equation} \theta = T_m \exp \left[-\frac{1}{2}\int _{R_c}^{R_p} \frac{\alpha g}{c_p} dr\right]. \end{equation} \tag{ 44 }$$

We adopt the relationship between Nusselt number, Rayleigh number, and viscosity contrast N across the conducting lid (Frank–Kamenetskii parameter) given by Solomatov & Moresi (2000);

$$\begin{equation} {\rm Nu} = a N^{-(1+\beta)}{\rm Ra}^\beta, \end{equation} \tag{ 45 }$$

where a ∼ 1, and boundary layer stability analysis shows that for a material with stress exponent n, β = n/(n + 2) (Solomatov 1995). We use β = 0.6.

Our numerical experiments revealed that the predicted history of the dynamo in a planet with an Earth-like interior is sensitive to (1) mantle viscosity, (2) the thermal expansivity of iron (related to its compressibility), (3) the thermal conductivity of Fe, (4) initial core temperature, and (5) the heat flow across the CMB. The first is constrained by the requirement to reproduce the size of and pressures in the Earth's core (Section 3.1). We parameterize the last by using different values of the critical Rayleigh number Ra_*. Although Ra_* is on the order of 10³ (Nimmo & Stevenson 2000; Turcotte & Schubert 2002), we consider a range of values that encapsulate the complex uncertainties of the thermochemical boundary layer at the CMB.

Mantle viscosity. Our viscosity law requires a solidus and a reference viscosity. We fit a Lindemann melting law to the solidus data of Zerr et al. (1998) for the lower mantle, using the Mg–perovskite equation of state to relate P to ρ: we found a best-fit Grüneisen parameter of 1.45. The reference viscosity of a hydrated Earth-like mantle was adjusted to achieve a present potential temperature θ of 1750 K. This is the present temperature beneath mid-ocean ridge basalts (1623 K; Herzberg et al. 2010) to which 130 K associated with the olivine–spinel transition in the deep mantle has been added (latent heat of reaction = 160 kJ kg⁻¹; Ito et al. 1971). The derived reference viscosity, η_* = 6 × 10¹⁹ Pa s, is somewhat above inferred values for the terrestrial asthenosphere where partial melting may occur (Solomatov & Moresi 2000; Dixon et al. 2004; Fjeldskaar & Cathles 2007; Bills et al. 2007; James et al. 2009). The viscosity of a dehydrated mantle is expected to be significantly higher. The effect of water on viscosity depends on the dominant creep mechanism, and the partitioning of water among the mineral phases that are present. Hydration of olivine in the upper mantle could lower the viscosity by a factor of ∼500 (Hirth & Kohlstedt 1996), but the degree of hydration and concomitant decrease in viscosity could be much less in the lower perovskite mantle (Bolfan-Casanova et al. 2000). We selected a value based on the requirement that a simulation of Venus (0.815 M_⊕ and T_p = 785 K) predicts no dynamo at the present time. The slow rotation of Venus cannot be solely responsible for the absent dynamo (Zhang & Zhang 1995). If Venus has had a stagnant lid over its entire history, then the required viscosity enhancement for the Venusian dynamo to cease by 4.5 Gyr is a factor of 6. If its lid was primarily mobile, then the required viscosity enhancement is a factor of 200. We conservatively adopt a factor of 10, recognizing that higher values are possible.

Fe thermal conductivity. Thermal conduction competes with advection in removing heat from the core. Higher thermal conductivity will suppress convection, and thus the dynamo, and will stratify some or all of the liquid core. The thermal conductivity of iron is expected to depend on pressure, temperature, and the presence of light elements, but not on temperature (Berman 1976). The value often used in core models is 50–60 W m⁻¹ K (Labrosse 2003). Thermal conductivity can be determined from measurements of electrical conductivity using the Wiedemann–Franz–Lorenz (WFL) law. Assuming the WFL law applies at high pressure, electrical conductivity measurement at 140 GPa (Keeler & Mitchell 1969) and 200 GPa (Bi et al. 2002) imply k ∼ 100 W m⁻¹ K⁻¹, cf., Anderson (1998). In contrast, Stacey & Loper (2007) argue for a value of 28–29 W m⁻¹ K⁻¹ based on the expectation that scattering of electrons by 3d orbitals (which contributes to resistivity) will increase faster with pressure than the densities of 3d and 4s states (which contribute to conductance). The electrical conductivity measurements of Matassov (1977) show that light elements can significantly reduce thermal conductivity (Anderson 1998), e.g., a factor of 2 reduction for Fe₃Si (Manga & Jeanloz 1996). A low value of k is an alternative mechanism to radiogenic heat to avoid core stratification (Stacey & Loper 2007).

Initial core temperature. We presume that the cores of more massive planets will be initially hotter because of the conversion of a greater amount of gravitational potential energy into heat. A hypothetical planet with a small core that grows concurrently by the descent of iron from a cold surface will have a final temperature T₀ ≈ (3GM_p)(10R_pc_p) ∼ 22100(M/M_⊕)^2/3 K. These core temperatures are unrealistic because much of the gravitational potential energy will be dissipated in and carried back to the surface by a partially molten silicate mantle. Instead, we set the initial temperature at the top of the core relative to the melting temperature of Mg–perovskite (MgSiO₃), effectively the liquidus temperature of the mantle, at the CMB. Core temperatures reaching this would result in the formation of a massive magma ocean whose low viscosity would quickly cool a hotter core. Alternatively, our initial state can be considered the epoch when any such magma ocean solidifies. The experimental data and ab initio calculations of the MgSiO₃ melting temperature (Stixrude et al. 2009) show that τ is well described to 136 GPa by the relationship τ≈ 5400(P/136 Gpa)^0.367. This is equivalent to a Lindemann law and γ = 1.44, cf. Akaogi & Ito (1993), and we simply extrapolate this to higher pressures. The mantle is given its steady-state temperature as its initial temperature.

Figure 3 plots the predicted inner core size and the surface magnetic field of the Earth at 4560 Myr for different values of the initial core temperature (from 100 K above to 500 K below the MgSiO₃ melting temperature), critical Rayleigh numbers between 460 and 3610, and three different values of the thermal conductivity (30, 35, and 40 W m⁻¹ K⁻¹). We calculate a present surface-average terrestrial field strength using a dipole moment of 8 × 10²² A m². To compare with the partially non-dipolar field of the Earth, we multiply our model predictions by 80% (Lowes 2007). Because 80% is the measured dipole fraction at the surface, and the fraction at the core will be smaller, our predictions should overestimate the strength of the Earth's current field. At k > 45 W m⁻¹ K⁻¹, the core becomes partly stratified, inconsistent with observations of magnetic field variability (Gubbins 2007; Stanley & Mohammadi 2008). At lower values of k, the present core size and magnetic field strength can be reproduced by combinations of a range of critical CMB Rayleigh numbers and initial core temperatures. We adopt k = 35 m⁻¹ K ⁻¹, Ra_* = 1100, and an initial core temperature of 275 K below the MgSiO₃ melting temperature at CMB pressure. We use these parameter values for all subsequent calculations.

Zoom In Zoom Out Reset image size

Figure 4 shows the calculated 4.56 Gyr evolution of a 1 M_⊕ planet. The mantle potential temperature decreases by 200 K, consistent with paleothermometry based on the chemistry of Archean and Proterozoic non-arc basalts (Herzberg et al. 2010). The temperature contrast across the CMB decreases until the onset of inner core formation, then rises slightly, as the core is heated by latent heat from the inner core, to a present value of around 1600 K, only slightly higher than current estimates (Tateno et al. 2009). Calculated mantle heat flow, radiogenic heat flow, and core heat flow are all consistent with current estimates (Buffett 2002; Lay et al. 2008), but the Urey number at 4.56 Gyr is ∼0.8, significantly higher than in some models (Korenaga 2008). The predicted evolution of the surface field (or the equivalent dipole moment) is comparable to that of other models (Labrosse 2007; Aubert et al. 2009; Breuer et al. 2009). It also does not conflict with paleomagnetic data, represented by a filtered sample of the IAGA paleointensity database (Biggin et al. 2009) as well as three recent measurements in 3.45 Ga rocks (Herrero-Bervera & Mojzsis 2009; Tarduno et al. 2010).

Zoom In Zoom Out Reset image size

Whether a core grows outward from the center of the planet, or inward from the top of the core depends on Δ, the ratio of the pressure–temperature slopes of the adiabat and solidus. This value will vary with depth in the core, but is prudently evaluated at the ICB or, when no core is present, at the center of the planet. For an adiabatic temperature profile,

$$\begin{equation} \Delta = \frac{2\left(3\gamma - 1\right)K_1 \left(1 + \ln (\rho _i/\rho _0)\right)}{\gamma \left[dP/d(\rho /\rho _1)\right]}. \end{equation} \tag{ 46 }$$

The zero-pressure compressibility K₁ appears as a linear factor in the EOS, therefore Δ ultimately depends only on the Grüneisen parameter γ, the pressure-derivative terms in the EOS, and the pressure at the ICB. Because gravity is comparatively low in the core, the pressure is relatively insensitive to the size of the core, and the total mass of the planet is the most important parameter. Δ will decrease with mass due to the increasing incompressibility of Fe (including the TFD effect) at high pressures. Δ < 1 for planets larger than 2.5–3 M_⊕ (Figure 5).

Zoom In Zoom Out Reset image size

Figure 6 shows the logarithmic ratio of the temperature to the melting point versus depth in the cores of 1–10 M_⊕ planets. An unimportant constant has been ignored in each case. The minima in these curves are where iron will first solidify. Above 2 M_⊕, the minimum shifts from the center (corresponding to inner core growth) to a position near or at the CMB. In these cases, iron "snow" will freeze out at or near the top of the core. Except under an exceptionally narrow range of conditions when the temperature coincides with the solidus throughout much of the core (Figure 6), iron crystals or slabs will first form at or near the top of the core and sink, and the resulting release of light elements and transport of latent heat will stabilize the liquid core and shut down convection (Williams 2009). This effect has been discussed for several small bodies in the solar system (Hauck et al. 2006; Stewart et al. 2007; Chen et al. 2008; Bland et al. 2008).

Zoom In Zoom Out Reset image size

Because there is no net release of latent heat, the core will continue to cool and iron will crystallize over an expanding depth range (Figure 6). Double diffusive convection can drive a "fingering" instability in a stratified fluid where two gradients, e.g., a thermal gradient and compositional gradient, oppose each other, but are associated with substantially different diffusivities (Huppert & Turner 1981; Huppert & Sparks 1984). Double diffusive convection has been included in a recent model of Mercury's dynamo (Manglik et al. 2010). In our case, the release and absorption of latent heat by crystallizing, sinking, and melting grains of iron will create a slightly subadiabatic temperature gradient, and buoyancy forces due to temperature will not play a role. Double diffusive convection will only occur if (1) the diffusion of the light elements is faster than the iron crystals and (2) the crystals are well-coupled to the fluid, i.e., there is a high density of nucleation sites and crystals are small. Even if convection ensues, it will do so in discrete layers each of which has a Rayleigh number of order the critical value (Huppert & Sparks 1984). Boundary layer theory for finite-amplitude convection (Turcotte & Schubert 2002) shows that the magnetic Reynolds number of each layer will be of order (κ/λ)Ra^2/3, or ∼10⁻³, well below that required to sustain a dynamo. If the crystals are not coupled to the fluid, no convective instability will occur and instead crystals will steadily settle through the fluid. It is not clear that a dynamo occurs under these conditions either.

For a given total mass, the mass of a planet's core may differ as a result of variation in the ratio of Fe to Si in the primordial disk (Neves et al. 2009), the efficiency of core formation (Elkins-Tanton & Seager 2008) or impact stripping of the mantle (Benz et al. 1988, 2007). The ability of impacts to remove substantial amounts of the mantle from Earth-mass and larger planets might be limited, and thus massive "Mercurys" may not be plausible (Marcus et al. 2010). Figure 7(a) shows the predicted dynamo history of Earth-mass planets with core sizes between 0.5 and 1.5 times that of the Earth. The average magnetic field is higher on planets with larger cores, mostly due to the smaller r⁻³ attenuation of the dipole field outside the core.

Zoom In Zoom Out Reset image size

The initial abundance of heat-producing radioactive nuclides (⁴⁰K, ²³²Th, ²³⁵U, and ²³⁸U) may also vary between planets and between systems as a result of chemical and thermal processing of solids in planet-forming disks, and galactic chemical evolution (Kite et al. 2009). K is a volatile element and the abundance of ⁴⁰K will vary with the location of a planet's accretion zone, with those planets (e.g., Mars) forming further out having a higher abundance (Lodders & Fegley 1997). The Galactic abundances of these isotopes relative to major planet-forming elements such as Si are predicted to have decreased substantially with time, and planets of a given age may accrete with substantially different inventories, but the difference between planets of different ages observed at the present time is much smaller (Kite et al. 2009). Figure 7(b) shows the predicted dynamo evolution for planets with 50% and 150% of the Earth's estimated initial inventory of radionuclides. The predicted surface field is slightly lower with a higher abundance of radionuclides, but is relatively insensitive to this parameter.

Figure 8 plots the predicted dynamo history of 1–10 M_⊕ Earth-like planets with and without plate tectonics, and surface temperatures of 288 K (Earth-like) and 1500 K (the solidus of basalt). Planets much larger than 1.5 M_⊕ with stagnant lids do not produce dynamos and even more massive planets have lower mantle temperatures that exceed the estimated peridotite solidus. Correct modeling of planets with internal magma oceans is beyond the scope of this study and we terminate those runs. Only planets smaller than 2.5 M_⊕ form inner cores (<2 M_⊕ for those with stagnant lids), and the elapsed time before an inner core appears, corresponding to a rapid increase in the strength of the magnetic field, increases with planet mass. The strength of the magnetic field at any given age decreases with planet mass. On Earth-mass planets with elevated surface temperatures, the strength of the magnetic field is higher and the interval before the onset of core formation is briefer.

Zoom In Zoom Out Reset image size

Planets more massive than 2 M_⊕ (like CoRoT-7b) have long-lived dynamos produced by core cooling, but only if they have a mobile lid. The duration of the dynamo decreases with planet mass and increases with surface temperature. The bottom right-hand panel of Figure 8 extends these trends out to a mass of 10 M_⊕.

For an inner core to form, enough sensible heat must be removed for the center temperature to reach the solidus. We obtain an approximate timescale t_c for this by setting Q_s = −Q_K and performing the appropriate integration. This gives

$$\begin{equation} t_c \approx \frac{M_c\Delta S}{8\pi kR_c} \left(1 - \frac{D^2}{R_c^2}\ln \frac{\tau _{0}({\rm Fe})}{\tau _{c}({\rm pv})}\right). \end{equation} \tag{ 47 }$$

where the solidi for Fe and perovskite are at central and CMB conditions, respectively. For k = 40 W m⁻¹ K⁻¹ and R_c ∼ M^1/4_p (Figure 2), the first factor is ≈2M^3/4_p Gyr. The second factor is of order unity: Both solidi increase with planet mass and so we expect t_c to increase slightly faster with planet mass than M^3/4_p. This is born out by our detailed calculations, i.e., a 2 M_⊕ planet takes more than twice as long to form an inner core than an Earth-mass planet (Figure 8).

A maximum time t_s for the complete solidification of the entire core, after the inner core starts growing, is obtained by dividing the total latent heat by the heat flow. (In our cases, the heat flow through the mantle tends to be close to or not much larger than the conduction through the core itself). In these cases, t_s is dependent on mass only through the pressure dependences of certain material properties of the core:

$$\begin{equation} t_s \approx \frac{c_p \Delta S}{2 k \alpha G}. \end{equation} \tag{ 48 }$$

For standard values, t_s ∼ 50 Gyr. The iron cores of rocky planets with Earth-like composition will never completely solidify, and solid inner cores never include more than about half of the total core mass in 10 Gyr. Thus, planets that do form inner cores (⩽2 M_⊕) will have dynamos that last longer than the main-sequence life of the parent star.

In planets with masses of 2–3 M_⊕, inner cores do not grow but instead iron condenses as "snow" at the CMB, shutting down convection and the dynamo within 10 Gyr. In yet more massive planets, the core never cools sufficiently for iron to condense. These planets have persistent dynamos produced by core cooling. However, the temperature contrast between the core and the mantle decays, and the heat flow across the CMB eventually falls to the core conduction value. At that point, ϕ = 0 and the dynamo ceases.

At a given planet mass, the surface field B_p increases with increasing core size, and a solid inner core formation forms earlier (Figure 7(a)). The former is a consequence of the cubic distance attenuation law and a shallower CMB, as well as increased heat flow across a thinner mantle. The latter is a result of more rapid cooling, but also higher pressure and a higher solidus at the planet's center. Dynamo history is much less sensitive to the initial abundance of long-lived radionuclides in the planet's mantle (Figure 7(b)) because of the self-regulating nature of viscosity-dependent mantle convection (Tozer 1972).

B_p generally decreases with planet mass. Equations (2) and (5) show why: the intensity of the field in the core B_c depends on the available entropy production per unit mass ϕ and only weakly on core size. All else being equal, ϕ will decrease with mass because the core of a larger planet has less surface area per unit mass and a thicker mantle through which to lose heat. If ϕ ∼ R⁻¹_c, this makes B_c insensitive to planet mass. The intensity of the field at the surface will decrease as R⁻³_c and thus will roughly scale as M^−3/4_p.

Higher surface temperature leads to elevated mantle temperatures and lower viscosities. Although the former lowers the temperature contrast across the CMB, the latter more than compensates because of the exponential decrease of viscosity with temperature. The outcome, increased magnetic field strength and dynamo lifetime, is clearly seen in Figure 8.

Our model predicts that CoRoT-7b has a dynamo at its estimated age of 1.3–2.2 Gyr (Léger et al. 2009), but only if it has a mobile lid. In this case, the dynamo is driven entirely by core cooling, and no inner core is present. In fact, it is the absence of iron solidification (as "snow") that allows a dynamo to persist. Whether or not plate tectonics operates on CoRoT-7b is debated (Valencia et al. 2007; O'Neill & Lenardic 2007; Kite et al. 2009), and the presence of a dynamo may be one means of testing this. However, a core dynamo is not the only mechanism that can produce a magnetic field on such an object (see Section 5.4).

High-pressure material properties. The properties of silicates and Fe at pressures of hundreds or thousands of GPa are very uncertain. Although EOS become dominated by electron degeneracy pressure at pressures well above 1 TPa, other parameters, namely k and α, are very important in determining whether a dynamo operates. Aitta (2006) has proposed that iron has a tricritical point at around 800 GPa at which separate solid and liquid phases disappear and the latent heat of fusion vanishes. These conditions are reached at the center of a 2.5 M_⊕ planet and at the CMB of a 6 M_⊕ planet (Figure 2).

Mantle structure. We assume a fully convecting mantle: a layered mantle will support a higher temperature contrast for the same heat flow (McNamara & van Keken 2000), and this would tend to suppress dynamo operation. The lowest 150 km of the Earth's mantle (within the "D" layer) may be in the stability field of a post-perovskite (ppv) phase (Oganov & Ono 2004) and ppv may dominate in sufficiently cool "super-Earth" mantles (Valencia et al. 2009). If a magma ocean existed at the base of the Hadean terrestrial mantle, release of latent heat during crystallization would have suppressed the cooling of Earth's core (Labrosse et al. 2007). Our calculations indicate that the mantles of planets with M ⩾ 2.5 M_⊕ and stagnant lids are partially molten. We do not attempt to model the thermal evolution of such objects but they are (at least initially) unlikely to support a dynamo.

Other effects of light element in the core. Our model does not include melting-point depression in the core by light elements, particularly sulfur (S). The effect of these is to delay inner core formation relative to the pure Fe case. In that sense, our calculations are optimistic in predicting dynamos. (This is in contrast to the case of small bodies where S can maintain a dynamo by delaying the freezing of the entire core.) However, we choose our initial core temperature based on the present size of the Earth's inner core and thus, relative to Earth, this calibration may cancel some of the effect of S. We include the production of buoyancy flux resulting from the exclusion of light elements from a growing solid core. However, we assume that light elements are uniformly mixed into the outer liquid core and do not include possible variations in light element concentrations with planet mass. The liquid silicate–liquid metal partition coefficient of sulfur increases with pressure (Li & Fei 2003), but no experimental data is yet available for pressures exceeding those at Earth's CMB. If the trend continues to higher mass (and for other elements), the power available for dynamos in the cores of super-Earths may be smaller than our model predicts.

Non-radiogenic sources of mantle heat. Dissipation of tidal strain could be an additional source of heat in the mantles of non-synchronously rotating planets close to their parent stars (Barnes et al. 2010). The timescale for rotational synchronization is a few megayears (Guillot et al. 1996) while that of orbital circularization can be gigayears (Rasio et al. 1996), but both are strongly dependent on semimajor axis. Sufficiently intense heating during this period could lead to a long-lived thermal anomaly in a mantle that influences the transport of heat much like a stagnant lid. We also ignore the potential complication of heating by late giant impacts (Roberts et al. 2009; Arkani-Hamed & Olson 2010).

Effect of the stellar magnetic field. In the absence of an ionosphere, the ambient magnetic field of the host star B_* could reach the planet's core: a sufficiently strong applied field will inhibit convective motion in transverse directions. This will occur when the Stuart number

$$\begin{equation} {\rm St} = \frac{\sigma B_*^2 R}{\rho v} \gg 1. \end{equation} \tag{ 49 }$$

Using the definition of Re_m, we rewrite this as

$$\begin{equation} \frac{\sigma B_*^2 R^2}{\rho \lambda } \gg {\mbox{Re}_m}. \end{equation} \tag{ 50 }$$

For typical values σ ∼ 5 × 10⁵ S m⁻¹, ρ ∼ 10⁴ kg m, λ ∼ 2 m² s⁻¹, L ∼ 3 × 10⁶ m, the condition becomes B_* ≫ 60 Re^1/2_m nT. For Re_m ∼ 10³, the requirement becomes B_* ≫ 2 μT. B_*∼ 100 μT near the surface of the Sun: assuming that the stellar magnetic field decreases with the square of the distance (because of interplanetary plasma), St ≫ 1 only within 0.1 AU. Thus, inhibition of convection could occur in the cores of the closest planets, although the increased temperature and conductivity of surface rocks at such proximity to the star may block the stellar field (see Section 5.4).

If a planet is embedded in an ambient magnetic field from its star that varies with orbital period P, eddy currents will generate an induced field if a layer is sufficiently thick (d) and conducting. The condition for a substantial induced field is (Stevenson 2003)

$$\begin{equation} d\sigma \gg \frac{1}{2\pi \mu _0 R_p}. \end{equation} \tag{ 51 }$$

For an Earth-size planet where d is in km and P in days, the required conductivity σ is ≫2Pd⁻¹(R_p/R_⊕)⁻¹ S m⁻¹. The conductivity of anhydrous molten silicates is 0.1–1 S m⁻¹ while that of molten carbonates is 3 orders of magnitude higher (Gaillard et al. 2008). A planet with a magma ocean only a few kilometers deep on a planet with a very short-period orbit such as CoRoT-7b (0.85 d) could have an induced field.

Uranus and Neptune have magnetic fields that are thought to be generated in a layer of high-pressure water or water ice (Cavazzoni et al. 1999). The ionic conductivity of planetary ices is about 2 × 10³ S m⁻¹ at pressures above 40 Gpa, conditions that are reached at 0.7–0.8 radii in these planets (Stanley & Bloxham 2006). We added liquid water and high-pressure ice (Ice VII) layers to the interior model described in Section 3.1 to calculate the pressure at the base of an upper water/ice mantle as a function of total planet mass and water mass fraction (Figure 9). A basal pressure of 40 Gpa is reached on planets of a few M_⊕ and ice fractions of ∼10%. Such an envelope of water could be easily accommodated by the observed radius (2.7 R_⊕) of GJ 1214b (Charbonneau et al. 2009) but not by the mass and radius of CoRoT-7b (Valencia et al. 2010).

Zoom In Zoom Out Reset image size

Our model of a planet consisting of a silicate mantle and pure Fe core reproduces the gross interior structure and thermal history of the Earth. Using the magnetic field strength scaling of Aubert et al. (2009), the predicted evolution of the field strength is also consistent with paleomagnetic inferences and the present value. When applied to Venus, with a stagnant lid, an elevated surface temperature, and enhanced mantle viscosity, our model correctly predicts the absence of a dynamo at present. We are thus confident in extrapolating the model to rocky planets with approximately Earth-like composition but with different mass, surface temperature, and/or mode of convection.

The model predicts that 1–2 M_⊕ planets with plate tectonics will sustain a dynamo for up to 7 Gyr by core cooling before an inner core forms and the dynamo intensifies at a time that depends on surface temperature, mass, and the relative size of the core. Planets in this mass range with stagnant lids will have dynamos only once an inner core begins to grow. Planets with higher surface temperatures can more easily sustain a dynamo because the resulting higher mantle temperatures lead to a lower viscosity and hence more rapid mantle convection and core cooling.

Planets larger than about 2 M_⊕ do not develop solid inner cores because the adiabat is steeper than the solidus. In such planets, iron freezes first at or near the CMB. The latent heat and any light elements released at the top of the core will not be able to drive a long-lived dynamo. In planets between 2–3 M_⊕, this leads to the premature shutdown of any dynamo within 10 Gyr.

In planets with masses >3 M_⊕, iron never begins to freeze in 10 Gyr. If these planets have plate tectonics, our model predicts that they have relatively weak dynamos driven solely by core cooling but which can last several gigayears, depending on mass and surface temperature. The mantles of planets with stagnant lids in this mass range are hot enough to prevent core cooling and dynamo operation, but they may also be partially molten. This last scenario is not addressed by our model. Our model predicts that CoRoT-7b may have a dynamo, but only if the planet has a mobile lid, and the strength of the field at its surface will be less than that of the present Earth (∼50 μT). Advances in our knowledge of the properties of planetary materials at very high pressure will improve our ability to make more quantitative predictions about the intensity and longevity of core dynamos in massive rocky exoplanets.