THE INFLUENCE OF THERMAL EVOLUTION IN THE MAGNETIC PROTECTION OF TERRESTRIAL PLANETS

Our one-dimensional model for the interior assumes a planet made by two well-differentiated chemically and mineralogically homogeneous shells: a rocky mantle made out of olivine and perovskite and a core made by iron plus other light elements.

The mechanical conditions inside the planet (pressure P, density ρ, and gravitational field g) are computed by simultaneously solving the continuity, Adams–Williamson, and hydrostatic equilibrium equations (Equations (1)–(4) in VAL06). For all planetary masses, we assume boundary conditions, ρ(r = R_p) = 4000 kg m⁻³ and P(r = R_p) = 0 Pa. For each planetary mass and core mass fraction CMF =M_core/M_p, we use an RK4 integrator and a shooting method to consistently compute the core R_c and planetary radius R_p.

For the sake of simplicity, we do not include in the interior model the two- or even three-layered structure of the mantle. Instead, we assume a mantle completely made of perovskite–postperovskite (ppv). This is the reason why we take ρ(r = R_p) = 4000 kg m⁻³ instead of the more realistic value of 3000 kg m⁻³. With a single layer and a realistic surface density, our model is able to reproduce the present interior properties of the Earth.

In all cases, we use the Vinet equation of state instead of the commonly used third-order Birch–Murnaghan equation (BM3). It is well known that the BM3 follows from a finite strain expansion and does not accurately predict the properties of the material for the typical pressures found in SEs, i.e., 100–1000 GPa for M_p = 1–10 M_⊕ (VAL06; Tachinami et al. 2011). We have ignored thermal corrections to the adiabatic compressibility, i.e., K_S(ρ, T) ≈ K_S(ρ, 300 K) + ΔK_s(T) (VAL06). This assumption allows us to decouple at runtime the CPU intensive calculation of the thermal profile from the mechanical structure at each time step in the thermal evolution integration.

Although we have ignored the "first-order" effect of the temperature in the mechanical structure, we have taken into account "second-order" effects produced by phase transitions inside the mantle and core. Using the initial temperature profile inside the mantle, we calculate the radius of transition from olivine to perovskite (neglecting the effect of an intermediate layer of wadsleyite). For that purpose, we use the reduced pressure function Π (Christensen 1985; Weinstein 1992; Valencia et al. 2007). For simplicity, the position of the transition layer is assumed to be constant during the whole thermal evolution of the planet. We have verified that this assumption does not significantly change the mechanical properties inside the planet, at least not at a level affecting the thermal evolution itself.

Inside the core, we continuously update the radius of the transition from solid to liquid iron (see below). For that purpose, we use the thermal profile computed at the previous time step. To avoid a continuous update of the mechanical structure, we assume in all the cases that during the transition from solid to liquid the density of iron changes by a constant factor Δρ = (ρ_s − ρ_l)/ρ_l. Here, the reference density of the solid ρ_s (when applied) is computed using the Vinet equation evaluated at a point in the very center of the planet.

Table 1 enumerates the relevant physical parameters of the interior structure and thermal evolution model.

Table 1. RTEM Parameters

Parameter	Definition	Value	Ref.
Bulk
CMF	Core mass fraction	0.325	⋅⋅⋅
Ts	Surface temperature	290 K	⋅⋅⋅
Ps	Surface pressure	0 bar	⋅⋅⋅
Inner core
⋅⋅⋅	Material	Fe	A
ρ0, K0, $K_0^{\prime }$, γ0, q, θ0	Equation of state parameters	8300 kg m−3, 160.2 GPa, 5.82, 1.36, 0.91, 998 K	A
kc	Thermal conductivity	40 W m−1 K−1	B
ΔS	Entropy of fusion	118 J kg−1 K−1	C
Outer core
⋅⋅⋅	Material	Fe(0, 8)FeS(0, 2)	A
ρ0, K0, $K_0^{\prime }$, γ0, q, θ0	Equation of state parameters	7171 kg m−3, 150.2 GPa, 5.675, 1.36, 0.91, 998 K	A
α	Thermal expansivity	1.4 × 10−6 K−1	D
cp	Specific heat	850 J kg−1 K−1	C
kc	Thermal conductivity	40 W m−1 K−1	B
κc	Thermal diffusivity	6.5 × 10−6 m2 s−1	E
ΔS	Entropy of fusion	118 J kg−1 K−1	C
adb	Adiabatic factor for Tc(t = 0)	0.7	⋅⋅⋅
ξc	Weight of Tc in core viscosity	0.4	⋅⋅⋅
Lower mantle
⋅⋅⋅	Material	pv+fmw	A
ρ0, K0, $K_0^{\prime }$, γ0, q, θ0	Equation of state parameters	4152 kg m−3, 223.6 GPa, 4.274, 1.48, 1.4, 1070 K	A
d, m, A, b, D0, mmol	Viscosity parameters	1 × 10−3 m, 2, 13.3, 12.33, 2.7 × 10−10 m2 s−1, 0.10039 kg mol−1	F
α	Thermal expansivity	2.4 × 10−6 K−1	D
cp	Specific heat	1250 J kg−1 K−1	C
km	Thermal conductivity	6 W m−1 K−1	C
κm	Thermal diffusivity	7.5 × 10−7 m2 s−1	E
ΔS	Entropy of fusion	130 J kg−1 K−1	C
Upper mantle
⋅⋅⋅	Material	Olivine	A
ρ0, K0, $K_0^{\prime }$, γ0, q, θ0	Equation of state parameters	3347 kg m−3, 126.8 GPa, 4.274, 0.99, 2.1, 809 K	A
B, n, E*, $\dot{\epsilon }$	Viscosity parameters	3.5 × 10−15 Pa−n s−1, 3, 430 × 103 j mol−1, 1 × 10−15 s−1	D
V*	Activation volume	2.5 × 10−6 m3 mol−1	F
α	Thermal expansivity	3.6 × 10−6 K−1	D
cp	Specific heat	1250 J kg−1 K−1	C
km	Thermal conductivity	6 W m−1 K−1	C
κm	Thermal diffusivity	7.5 × 10−7 m2 s−1	E
θ	Potential temperature	1700 K	F
χr	Radioactive heat correction	1.253	⋅⋅⋅

Note. Sources: (A) VAL06, (B) Nimmo (2009), (C) Gaidos et al. (2010), (D) Tachinami et al. (2011), (E) Ricard (2009), (F) Stamenković et al. (2011).

Download table as:  ASCIITypeset image

In order to compute the thermal evolution of the core, we solve the equations of energy and entropy balance (Labrosse et al. 2001; Nimmo 2009):

$$\begin{equation} Q_c = Q_s + f_i(Q_g + Q_l) \end{equation} \tag{ 7 }$$

$$\begin{equation} \Phi = E_l + f_i(E_s + E_g) - E_k. \end{equation} \tag{ 8 }$$

Here, Q_c is the total heat flowing through the CMB and Φ is the total entropy dissipated in the core. E_s and Q_s are the entropy and heat released by the secular cooling, E_g and Q_g are the contribution to entropy and heat due to the redistribution of gravitational potential when light elements are released at the liquid–solid interface (buoyant energy), E_l and Q_l are the entropy and heat released by the phase transition (latent heat), and E_k is a term accounting for the sink of entropy due to the conduction of heat along the core. We have avoided the terms coming from radioactive and pressure heating because their contribution is negligible at the typical conditions inside SEs (Nimmo 2009). As long as the buoyant and latent entropy and heat are only present when a solid inner core exists, we have introduced a boolean variable f_i that turns on these terms when the condition for the solidification of the inner core arises.

The terms in the energy and entropy balance are a function of the time derivative of the temperature profile ∂T(r, t)/∂t (for detailed expressions of these terms, see Table 1 in Nimmo 2009). As an example, the secular heat and entropy are given by

$$\begin{eqnarray} Q_s & = & -\int \rho c_p\frac{\partial T(r,t)}{\partial t}dV, \nonumber \\ E_s & = & -\int \rho c_p \left[\frac{1}{T_c(t)} - \frac{1}{T(r,t)}\right]\frac{\partial T(r,t)}{\partial t}dV. \end{eqnarray} \tag{ 9 }$$

Here, c_p is the specific heat of the core alloy and T_c is the temperature at the CMB. If we assume that the temperature profile of the core does not change during the thermal evolution, we can write temperature as T(r, t) = f_c(r)T_c(t). Here, f_c(r) is the core temperature radial profile that we will assume adiabatic (see below). It should be noted again that T_c(t) = T(r = R_c, t).

With this assumption the energy balance equation (7) can be written as a first-order differential equation on T_c:

$$\begin{equation} Q_c=M_c [C_s + f_i(C_g + C_l)] \frac{ d T_c }{d t}, \end{equation} \tag{ 10 }$$

where M_c is the total mass of the core and C_s, C_g, and C_l are core bulk heat capacities which can be expressed as volumetric integrals of the radial profile f_c(r). In this equation, the total heat Q_c is intrinsically a function of T_c and should be computed independently (see below).

Using a simple exponential fit for the core density, as proposed by Labrosse et al. (2001), the adiabatic temperature profile can be approximated as (Labrosse 2003)

$$\begin{equation} f_c(r) = \exp \left(\frac{R_c^2 - r^2}{D_c^2}\right), \end{equation} \tag{ 11 }$$

where $D_c=\sqrt{3 c_p/2\pi \alpha \rho _c G}$ is the scale height of temperature, α is the isothermal expansivity (assumed for simplicity constant along the core), and ρ_c is the density at core center. Using this fit, the bulk secular heat capacity C_s ≡ Q_s/(M_c dT_c/dt) can be obtained from Equation (9):

$$\begin{equation} C_s=-4\pi \int _0^{R_c}\rho (r)c_p\exp \left(\frac{R_c^2 - r^2}{D_c^2}\right) r^2 dr. \end{equation} \tag{ 12 }$$

Analogous expressions for C_g and C_l are obtained from the definition of Q_g and Q_l as given in Table 1 of Nimmo (2009).

The total heat released by the core Q_c(T_c) is calculated here using the boundary layer theory (BLT; see, e.g., Stevenson 1983). Under this approximation Q_c is given by (Ricard 2009)

$$\begin{equation} Q_c = {4\pi}\!{R}_c k_m \Delta T_{\rm CMB} {\rm Nu}_c, \end{equation} \tag{ 13 }$$

where k_m is the thermal conductivity of the lower mantle, ΔT_CMB = T_c − T_l is the temperature contrast across the CMB, T_l is the lower mantle temperature, and Nu_c ≈ (Ra_c/Ra_*)^1/3 is the Nusselt number at the core (Schubert et al. 2001). The critical Rayleigh number Ra_* is a free parameter in our model (see Table 1).

The local Rayleigh number Ra_c at the CMB is calculated under the assumption of a boundary heated from below (Ricard 2009),

$$\begin{equation} {\rm Ra}_c = \frac{\rho \; g\; \alpha \;\Delta T_{\rm CMB}(R_p - R_c)^3}{\kappa _c \eta _c}, \end{equation} \tag{ 14 }$$

where g is the gravitational field and κ_c the thermal diffusivity at the CMB. The value of the dynamic viscosity η_c, which is strongly dependent on temperature, could be suitably computed using the so-called film temperature (see Manga et al. 2001 and references therein). This temperature could be computed in general as a weighted average of the temperatures at the boundaries,

$$\begin{equation} T_{{\eta }c}=\xi _cT_c + (1-\xi _c) T_l, \end{equation} \tag{ 15 }$$

where the weighting coefficient ξ_c is a free parameter whose value is chosen in order to reproduce the thermal properties of the Earth (see Table 1).

To model the formation and evolution of the solid inner core, we need to compare at each time the temperature profile with the iron solidus. We use here the Lindemann law as parameterized by VAL06:

$$\begin{equation} \frac{ \partial \log \tau }{\partial \log \rho } = 2 \left[ \gamma - \delta (\rho) \right]. \end{equation} \tag{ 16 }$$

Here, δ(ρ) ≈ 1/3 and γ is an effective Grüneisen parameter that, for simplicity, is assumed to be constant. To integrate this equation, we use the numerical density profile provided by the interior model and the reference values ρ₀ = 8300 kg m⁻³ (pure iron) and τ₀ = 1808 K.

The central temperature T(r = 0, t) and the solidus at that point τ(r = 0) are compared at each time step. When T(0, t_ic) ≈ τ(0) (t_ic is the time of inner core formation) we turn on the buoyant and latent heat terms in Equations (7) and (8), i.e., set f_i = 1, and continue the integration including these terms. The radius of the inner core at times t > t_ic is obtained by solving the equation proposed by Nimmo (2009) and further developed by Gaidos et al. (2010),

$$\begin{equation} \frac{{dR}_{\rm ic} }{d t} =-\frac{D_c^2}{2R_{\rm ic}(\Delta - 1)}\frac{1}{T_c} \frac{ d T_c }{d t}. \end{equation} \tag{ 17 }$$

Here, Δ is the ratio between the gradient of the solidus (Equation (16)) and the actual temperature gradient T_c(t)f_c(r) as measured at R_ic(t).

When the core cools down below a given level the outer layers start to stratificate. Here, we model the effect of stratification by correcting the radius and temperature of the core following the prescription by Gaidos et al. (2010). When stratified the effective radius of the core is reduced to R_⋆ (Equation (27) in Gaidos et al. 2010) and the temperature at the core surface is increase to T_⋆ (Equation (28) in Gaidos et al. 2010). The stratification of the core reduces the height of the convective shell which leads to a reduction of the Coriolis force potentially enhancing the intensity of the dynamo-generated magnetic field.

The estimation of the dynamo properties requires the computation of the available convective power Q_conv. Q_conv is calculated here assuming that most of the dissipation occurs at the top of the core. Under this assumption,

$$\begin{equation} Q_{\rm conv}(t)\approx \Phi (t) T_c(t), \end{equation} \tag{ 18 }$$

where the total entropy Φ is computed from the entropy balance (Equation (8)) using the solution for the temperature profile T_c(t)f_c(r).

When Φ(t) becomes negative, i.e., E_l + E_s + E_g < E_k in Equation (8), Q_conv also becomes negative and convection is no longer efficient to transport energy across the outer core. Under this condition, the dynamo is shut down. The integration stops when this condition is fulfilled at a time we label as the dynamo lifetime t_dyn.

One of the novel features of our thermal evolution model is that we treat mantle thermal evolution with a similar formalism as that described before for the metallic core.

The energy balance in the mantle can be written as

$$\begin{equation} Q_m = \chi _r Q_r + Q_s + Q_c. \end{equation} \tag{ 19 }$$

Here, Q_m is the total heat flowing out through the surface boundary (SB), Q_r is the heat produced in the decay of radioactive nuclides inside the mantle, Q_s is the secular heat, and Q_c is the heat coming from the core (Equation (13)).

We use here the standard expressions and parameters for the radioactive energy production as given by Kite et al. (2009). However, in order to correct for the non-homogeneous distribution of radioactive elements in the mantle, we introduce a multiplicative correction factor χ_r. Here, we adopt χ_r = 1.253 that fits the Earth properties well. We have verified that the thermal evolution is not too sensitive to χ_r and have assumed the same value for all planetary masses.

The secular heat in the mantle is computed using an expression analogous to Equation (9). As in the case of the core, we assume that the temperature radial profile does not change during the thermal evolution. Under this assumption, the temperature profile in the mantle can be also written as T(r, t) = T_m(t)f_m(r). In this case, T_m(t) = T(r = R_p, t) is the temperature just below the SB layer (see Figure 1).

Zoom In Zoom Out Reset image size

Assuming an adiabatical temperature profile in the mantle, we can also write:

$$\begin{eqnarray*} &&f_m(r) = \exp \left(\frac{R_p^2 - r^2}{D_m^2}\right), \end{eqnarray*}$$

in analogy to the core temperature profile (Equation (11)). In this expression, D_m is the temperature scale height for the mantle which is related to the density scale height L_m through $D_m^2=L_m^2/\gamma$ (Labrosse 2003). In our simplified model, we take the values of the density at the boundaries of the mantle and analytically obtain an estimate for L_m and hence for D_m.

The energy balance in the mantle is balanced when we independently calculate the heat Q_m at the SB as a function of T_m. In this case, whether or not mobile lids are present plays an important role in determining the efficiency with which the planet gets rid of the heat coming from the mantle. In the mobile lid regime, we assume that the outer layer is fully convective and use the BLT approximation to calculate Q_m,

$$\begin{equation} Q_{m}^{\rm ML} = \frac{{4\pi}{R}_p^2 k_m \Delta T_m {\rm Nu}_m}{R_p - R_c}, \end{equation} \tag{ 20 }$$

where ΔT_m = T_m − T_s is the temperature contrast across the SB and T_s is the surface temperature. Since we are studying the thermal evolution of habitable planets, in all cases we assume T_s = 290 K. Planetary interior structures and thermal evolution are not too sensitive to surface temperature. We have verified that the results are nearly the same for the surface temperature in the range of 250–370 K. In the mobile lid regime, Nu_m obeys the same relationship with the critical Rayleigh number as in the core. In this case, however, we compute the local Rayleigh number under the assumption of material heated from inside (Gaidos et al. 2010),

$$\begin{equation} {\rm Ra}_m = \frac{\alpha g \rho ^2 H (R_p - R_c)^5}{k_m \kappa _m \eta _m}, \end{equation} \tag{ 21 }$$

where H = (Q_r + Q_c)/M_m is the density of heat inside the mantle, k_m and κ_m are the thermal conductivity and diffusivity, respectively, and η_m is the upper mantle viscosity.

In the stagnant lid regime, the SB provides a rigid boundary for the heat flux. In this case, we adopt the approximation used by Nimmo & Stevenson (2000):

$$\begin{equation} Q_{m}^{\rm SL} = 4\pi R_p^2 \frac{k_m}{2} \left(\frac{\rho g \alpha }{\kappa _m \eta _m} \right) ^{1/3}\Gamma ^{-4/3}. \end{equation} \tag{ 22 }$$

Here, Γ ≡ −∂ln η_m/∂T_m measures the viscosity dependence on temperature evaluated at the average mantle pressure.

With all these elements at hand the energy balance at Equation (19) is finally transformed into an ordinary differential equation for the upper mantle temperature T_m(t),

$$\begin{equation} Q_m = \chi _r Q_r + Q_c + C_m \frac{ d T_m }{d t}, \end{equation} \tag{ 23 }$$

where C_m is the bulk heat capacity of the upper mantle which is calculated with an expression analogous to Equation (12).

In order to solve the coupled differential equations (10), (17), and (23), we need to choose a proper set of initial conditions.

The initial value of the upper mantle temperature is chosen using the prescription by Stamenković et al. (2011). According to this prescription, T_m(t = 0) is computed by integrating the pressure-dependent adiabatic equation up to the average pressure inside the mantle 〈P_m〉,

$$\begin{equation} T_{m}(t=0) = \theta \exp \left(\int _0^{\langle P_m\rangle }\frac{\gamma _0}{K_s(P^{\prime })}{dP}^{\prime }\right). \end{equation} \tag{ 24 }$$

Here, θ = 1700 K is a potential temperature which is assumed to be the same for all planetary masses (Stamenković et al. 2011). Using T_m(t = 0) and the adiabatic temperature profile, we can obtain the initial lower mantle temperature T_l(t = 0).

The initial value of the core temperature T_c(t = 0) is one of the most uncertain parameters in thermal evolution models. Although its actual value or its dependence on the formation history and planetary mass is not known, it is reasonable to start the integration of a simplified thermal evolution model when the core temperature is of the same order as the melting point for MgSiO₃ at the lower mantle pressure. A small arbitrary temperature contrast against this reference value (Gaidos et al. 2010; Tachinami et al. 2011) or more complicated mass-dependent assumptions (Papuc & Davies 2008) have been used in previous models to set the initial core temperature. We use here a simple prescription that agrees reasonably well with previous attempts and provides a unified expression that could be used consistently for all planetary masses.

According to our prescription, the temperature contrast across the CMB is assumed to be proportional to the temperature contrast across the whole mantle, i.e., ΔT_CMB = _adbΔT_adb = _adb(T_m − T_l). We have found that the thermal evolution properties of Earth are reproduced when we set _adb = 0.7.

Using this prescription, the initial core temperature is finally calculated using

$$\begin{equation} T_c(t=0)=T_l+\epsilon _{\rm adb}\Delta T_{\rm adb}. \end{equation} \tag{ 25 }$$

We have observed that the value of the T_c(t = 0) obtained with this prescription is very close to the perovskite melting temperature at the CMB for all the planetary masses studied here. This result shows that although our criterion is not particularly better physically rooted than those used in previous models, it still relies in just one free parameter, i.e., the ratio of mantle and CMB contrasts _adb.

One of the most controversial aspects and probably the largest source of uncertainties in thermal evolution models is the calculation of the rheological properties of silicates and iron at high pressures and temperatures. A detailed discussion on this important topic is beyond the scope of this paper. An up-to-date discussion and analysis of the dependence on pressure and temperature of viscosity in SEs and its influence in thermal evolution can be found in the recent works by Tachinami et al. (2011) and Stamenković et al. (2011).

We use here two different models to calculate viscosity under different ranges of temperatures and pressures. For the high pressures and temperatures of the lower mantle, we use a Nabarro–Herring model (Yamazaki & Karato 2001),

$$\begin{equation} \eta _{\rm NH}(P,T) = \frac{R_gd^m}{D_0 A m_{{\rm mol}}}T\rho (P,T)\exp \left(\frac{b\;T_{\rm melt}(P)}{T} \right). \end{equation} \tag{ 26 }$$

Here, R_g = 8.31 J mol⁻¹ K⁻¹ is the gas constant, d is the grain size, m is the growing exponent, A and b are free parameters, D₀ is the pre-exponential diffusion coefficient, m_mol is the molar density of perovskite, and T_melt(P) is the melting temperature of perovskite that can be computed with the empirical fit:

$$\begin{eqnarray*} &&T_{\rm melt}(P) = \sum _{i=0}^4 a_i \cdot P^i. \end{eqnarray*}$$

All the parameters used in the viscosity model, including the expansion coefficients a_i in the melting temperature formula, were taken from the recent work by Stamenković et al. (2011). The Nabarro–Herring formula allows us to compute η_c = η_NH(T_ηc), where T_ηc is the film temperature computed using the average in Equation (15).

The upper mantle has a completely different mineralogy and it is under the influence of lower pressures and temperatures. Although previous works have used the same model and parameters to calculate viscosity across the whole mantle (Stamenković et al. 2011, for example, use the perovskite viscosity parameters also in the olivine upper mantle), we have found here that using a different rheological model in the upper and lower mantles avoids under- and overestimation, respectively, of the value of viscosity that could have a significant effect on the thermal evolution.

In the upper mantle, we find that using an Arrhenius-type model leads to better estimates of viscosity than that obtained using the Nabarro–Herring model. In the upper mantle, the Nabarro–Herring formula (which is best suited to describe the dependence on viscosity at high pressures and temperatures) leads to huge underestimations of viscosity in that region. In the case of the Earth, this underestimation produces values of the total mantle too high compared to that observed in our planet making impossible to fit the thermal evolution of a simulated Earth.

For the Arrhenius-type formula, we use the same parameterization given by Tachinami et al. (2011):

$$\begin{equation} \eta _A(P,T) = \frac{1}{2} \left[ \frac{1}{B^{1/n}}\exp \left(\frac{E^* + PV^*}{nR_gT} \right) \right] \dot{\epsilon }^{(1-n)/n}, \end{equation} \tag{ 27 }$$

where $\dot{\epsilon }$ is the strain rate, n is the creep index, B is the Barger coefficient, and E* and V* are the activation energy and volume. The values assumed here for these parameters are the same as that given in Table 4 of Tachinami et al. (2011) except for the activation volume whose value we assume here V* = 2.5 × 10⁻⁶ m³ mol⁻¹. Using the formula in Equation (27), the upper mantle viscosity is computed as η_m = η_A(〈P_m〉, T_m).

A summary of the parameters used by our interior and thermal evolution models is presented in Table 1. The values listed in Column 3 define what we will call the reference thermal evolution model (RTEM). These reference values have been mostly obtained by fitting the present interior properties of the Earth and the global features of its thermal, dynamo, and magnetic field evolution (time of inner core formation and present values of R_ic, Q_m and surface magnetic field intensity). For the stagnant lid case, we use, as suggested by Gaidos et al. (2010), the values of the parameters that globally reproduce the present thermal and magnetic properties of Venus.

Figure 2 shows the results of applying the RTEM to a set of hypothetical TPs in the mass range M_p = 0.5–6 M_⊕.

Zoom In Zoom Out Reset image size

The surface temperature and hence "first-order" habitability of a planet depends on three basic factors: (1) the fundamental properties of the star (luminosity L_⋆, effective temperature T_⋆, and radius R_⋆), (2) the average star–planet distance (distance to the HZ), and (3) the commensurability of planetary rotation and orbital period (tidal locking). These properties should be properly modeled in order to assess the degree of star–planet interaction which is critical for determining the magnetic protection.

The basic properties of main-sequence stars of different masses and metallicities have been studied for decades and are becoming critical for assessing the actual properties of newly discovered exoplanets. The case of low-mass main-sequence stars (GKM) are particularly important for providing the properties of the stars with the highest potential to harbor habitable planets with evolved and diverse biospheres.

In this work, we will use the theoretical results of Baraffe et al. (1998, hereafter BAR98) that predict the evolution of different metallicities main-sequence GKM stars. We have chosen from that model those results corresponding to the case of solar metallicity stars. We have disregarded the fact that the basic stellar properties actually evolve during the critical period where magnetic protection will be evaluated, i.e., t = 0.5–3 Gyr. To be consistent with the purpose of estimating upper limits of magnetic protection, we took the stellar properties as provided by the model at the highest end of the time interval, i.e., t = 3 Gyr. Since luminosity increases with time in GKM stars, this assumption guarantees the largest distance of the HZ and hence the lowest effects of the stellar insolation and the stellar wind.

In order to estimate the HZ limits, we use the recently updated values calculated by Kopparapu et al. (2013). In particular, we use the interpolation formula in Equation (2) and coefficients in Table (2) to compute the most conservative limits of recent Venus and early Mars. The limits calculated for the stellar properties assumed here are depicted in Figure 5.

Zoom In Zoom Out Reset image size

The orbital and rotational properties of planets at close-in orbits are strongly affected by gravitational and tidal interactions with the host star. Tidal torque dampens the primordial rotation and axis tilt leaving the planet in a final resonant equilibrium where the period of rotation P becomes commensurable with the orbital period P_o,

$$\begin{equation} P:P_o=n:2. \end{equation} \tag{ 33 }$$

Here, n is an integer larger than or equal to 2. The value of n is determined by multiple dynamic factors, the most important being the orbital planetary eccentricity (Leconte et al. 2010; Ferraz-Mello et al. 2008; Heller et al. 2011). In the solar system, the tidal interaction between the Sun and Mercury has trapped the planet in a 3:2 resonance. In the case of GL 581d, detailed dynamic models predict a resonant 2:1 equilibrium state (Heller et al. 2011), i.e., the rotation period of the planet is a half of its orbital period.

Although in general estimating the time required for the "tidal erosion" is very hard given the large uncertainties in the key physical parameters involved (see Heller et al. 2011 for a detailed discussion), the maximum distance a_tid at which a solid planet in a circular orbit becomes tidally locked before a given time t can be roughly estimated by (Peale 1977)

$$\begin{equation} a_{\rm tid}(t) = 0.5\,{\rm AU}\,\left[\frac{({M_\star /M_\odot })^2 {P}_{\rm prim}} {{Q}}\right]^{1/6} t^{1/6}. \end{equation} \tag{ 34 }$$

Here, the primordial period of rotation P_prim should be expressed in hours, t in Gyr, and Q is the dimensionless dissipation function. For the purposes of this work, we assume a primordial period of rotation P_prim = 17 hr (Varga et al. 1998; Denis et al. 2011) and a dissipation function Q ≈ 100 (Henning et al. 2009; Heller et al. 2011).

In Figure 5, we summarize the properties of solar metallicity GKM main-sequence stars provided by the BAR98 model and the corresponding limits of the HZ and tidal locking maximum distance. The properties of the host stars of the potentially habitable SEs already discovered, GL 581d, GJ 667Cc, and HD 40307g, are also highlighted in this figure.

The stellar wind and CRs pose the highest risks for a magnetically unprotected potentially habitable TP. The dynamic pressure of the wind is able to obliterate an exposed atmosphere, especially during the early phase of stellar evolution (Lammer et al. 2003), and energetic stellar CRs could pose a serious risk to any form of surface life directly exposed to them (Grießmeier et al. 2005).

The last step in order to estimate the magnetospheric properties and hence the level of magnetic protection is predicting the stellar wind properties for different stellar masses and as a function of planetary distance and time.

There are two simple models used to describe the spatial structure and dynamics of the stellar wind: the pure hydrodynamical model originally developed by Parker (1958) that describes the wind as a non-magnetized, isothermal, and axially symmetric flux of particles (hereafter Parker's model), and the more detailed albeit simpler magneto-hydrodynamic model originally developed by Weber & Davis (1967) that takes into account the effects of stellar rotation and treats the wind as a magnetized plasma.

It has been shown that Parker's model reliably describes the properties of the stellar winds in the case of stars with periods of rotation of the same order as the present solar value, i.e., P ∼ 30 days (Preusse et al. 2005). However, for rapidly rotating stars, i.e., young stars and/or active dM stars, the isothermal model underestimates the stellar wind properties by almost a factor of two (Preusse et al. 2005). For the purposes of scaling the properties of the planetary magnetospheres (Equations (3)–(6)), an underestimation of the stellar wind dynamic pressure of that size will give us values of the key magnetospheric properties that will be off by 10%–40% of the values given by more detailed models. Magnetopause fields that have the largest uncertainties will be underestimated by ∼40%, while standoff distances and polar cap areas will be, respectively, under- and overestimated by just ∼10%.

According to Parker's model, the stellar wind average particle velocity v at distance d from the host star is obtained by solving Parker's wind equation (Parker 1958):

$$\begin{equation} u^2-\log u=4 \log \rho + \frac{4}{\rho } - 3, \end{equation} \tag{ 35 }$$

where u = v/v_c and ρ = d/d_c are the velocity and distance normalized with respect to $v_c=\sqrt{k_B T/m}$ and d_c = GM_⋆m/(4k_BT) which are, respectively, the local sound velocity and the critical distance where the stellar wind becomes subsonic. T is the temperature of the plasma which, in the isothermal case, is assumed to be constant at all distances and equal to the temperature of the stellar corona. T is the only free parameter controlling the velocity profile of the stellar wind.

The number density n(d) is calculated from the velocity using the continuity equation:

$$\begin{equation} n(d)=\frac{\dot{M_{\star }}}{4\pi d^2 v(d) m}. \end{equation} \tag{ 36 }$$

Here, $\dot{M_{\star }}$ is the stellar mass-loss rate, which is a free parameter in the model.

To calculate the evolution of the stellar wind, we need a way to estimate the evolution of the coronal temperature T and the mass-loss rate $\dot{M_{\star }}$.

Using observational estimates of the stellar mass-loss rate (Wood et al. 2002) and theoretical models for the evolution of the stellar wind velocity (Newkirk 1980), Grießmeier et al. (2004) and Lammer et al. (2004) developed semiempirical formulae to calculate the evolution of the long-term-averaged number density and velocity of the stellar wind for main-sequence stars at a given reference distance (1 AU):

$$\begin{equation} v_{1\,\rm {AU}}(t)=v_0 \left(1+\frac{t}{\tau }\right)^{\alpha _v} \end{equation} \tag{ 37 }$$

$$\begin{equation} n_{1\,\rm {AU}}(t)=n_0 \left(1+\frac{t}{\tau }\right)^{\alpha _n}. \end{equation} \tag{ 38 }$$

Here, α_v = −0.43, α_n = −1.86 ± 0.6, and τ = 25.6 Myr (Grießmeier et al. 2009). The parameters v₀ = 3971 km s⁻¹ and n₀ = 1.04 × 10¹¹ m⁻³ are estimated from the present long-term averages of the solar wind as measured at the distance of the Earth n(4.6 Gyr, 1 AU, 1 M_☉) = 6.59 × 10⁶ m⁻³ and v(4.6 Gyr, 1 AU, 1 M_☉) = 425 km s⁻¹ (Schwenn 1990).

Using these formulae, Grießmeier et al. (2007a) devised a clever way to estimate consistently T(t) and $\dot{M_{\star }}(t)$ in Parker's model and hence we are able predict the stellar wind properties as a function of d and t. For the sake of completeness, we summarize this procedure here. For further details, see Section 2.4 in Grießmeier et al. (2007a).

For a stellar mass M_⋆ and time t, the velocity of the stellar wind at d = 1 AU, v_1 AU, is calculated using Equation (37). Replacing this velocity in Parker's wind equation for d = 1 AU, we numerically find the temperature of the corona T(t). This parameter is enough to provide the whole velocity profile v(t, d, M_⋆) at time t. To compute the number density, we need the mass-loss rate for this particular star and at this time. Using the velocity and number density calculated from Equations (37) and (38), the mass-loss rate for the Sun $\dot{M_\odot }$ at time t and d = 1 AU can be obtained:

$$\begin{equation} \dot{M}_\odot (t)=4\pi (1\,{\rm AU})^2\;m\; n_{1\,\rm {AU}}({t}) v_{1\,\rm {AU}}({t}). \end{equation} \tag{ 39 }$$

Assuming that the mass-loss rate scales simply with the stellar surface area, i.e., $\dot{M}_\star (t)=\dot{M}_\odot (t)(R_\star /R_\odot)^2$, the value of $\dot{M_{\star }}$ can finally be estimated. Using v(t, d, M_⋆) and $\dot{M_{\star }}$ in the continuity equation (36), the number density of the stellar wind n(t, d, M_⋆) is finally obtained.

The value of the stellar wind dynamic pressure P_dyn(t, d, M_⋆) = m n(t, d, M_⋆) v(t, d, M_⋆)² inside the HZ of four different stars as computed using the procedure described previously is plotted in Figure 6.

Zoom In Zoom Out Reset image size

It is important to stress here that for stellar ages t ≲ 0.7 Gyr the semiempirical formulae in Equations (37) and (38) are not longer reliable (Grießmeier et al. 2007a). These equations are based on the empirical relationship observed between the X-ray surface flux and the mass-loss rate $\dot{M_{\star }}$ (Wood et al. 2002, 2005) which has been reliably obtained only for ages t ≳ 0.7 Gyr. However, Wood et al. (2005) have shown that an extrapolation of the empirical relationship to earlier times overestimates the mass-loss rate by a factor of 10–100. At times t ≲ 0.7 Gyr and over a given magnetic activity threshold, the stellar wind of main-sequence stars seems to be inhibited (Wood et al. 2005). Therefore, the limit placed by observations at t ≈ 0.7 Gyr is not simply an observational constraint but could mark the time where the early stellar wind also reaches a maximum (J. L. Linsky 2012, private communication). This fact suggests that at early times the effect of the stellar wind on the planetary magnetosphere is much lower than normally assumed. Hereafter, we will assume that intrinsic PMF is strong enough to protect the planet, at least until the maximum of the stellar wind is reached at t ≈ 0.7 Gyr, and focus on the stellar wind and magnetosphere properties for times larger than this.

Combining the model of magnetosphere evolution developed here with models of thermal and non-thermal atmospheric escape, it would be possible to estimate the mass-loss rate from atmospheres of magnetized and unmagnetized potentially habitable planets. This is a fundamental goal to be pursued in the near future if we want to assess the actual habitability of TPs in the HZ of their host stars presently known and those discovered in the future. The complex interaction between an inflated atmosphere and its protective magnetosphere and large uncertainties in the surface fluxes of atmospheric gasses that compensate the loss of volatiles induced by the action of the stellar wind render this goal hard to achieve at present. Despite these limitations, we can still make order of magnitude estimations based on our own results and the mass-loss rate computed, for example, in the recent works by Tian et al. (2008), Tian (2009), and Lammer et al. (2012).

Atmospheric thermal mass loss induced by the exposition to X-rays and EUV stellar radiation (XUV) have been estimated for the case of Earth-like N₂-rich atmospheres (Watson et al. 1981; Kulikov et al. 2006; Tian et al. 2008) and dry Venus-like CO₂-rich atmospheres (Tian 2009; Lammer et al. 2012). One critical property of an inflated atmosphere is essential to evaluate the exposition of such atmospheres to further non-thermal processes: the radius of the exobase R_exo. R_exo is defined as the distance where the mean-free path of atmospheric particles could be comparable to the size of the planet. When the radius of the exobase is comparable to or larger than the magnetic standoff distance R_S, we will say that the planet is unmagnetized. Under these conditions, the gases escaping from the exosphere will be picked up by the stellar wind and lost to the interplanetary space. On the other hand, if the exobase is well inside the magnetosphere (which is the case of the Earth today) atmospheric gasses escaping thermally from the exosphere could stay trapped by the magnetic field forming a plasmasphere. Planets under this condition will be magnetically protected and the mass-loss rate is expected to be much lower than for unmagnetized planets.

Using the conservative estimation of the X and EUV luminosities of main-sequence stars given by Garcés et al. (2011), we have estimated the XUV flux at the top of the atmospheres of GL 581d, GJ 667Cc, and HD 40307g during the first critical gigayear of planetary evolution. The planet that received the minimum amount of XUV radiation is HD 40307g with F_XUV = 10–35 PEL (1 PEL = 0.64 erg cm⁻² s⁻¹ is the present Earth value; Judge et al. 2003; Guinan et al. 2009). GL 581d was exposed in the first gigayear to a flux of F_XUV = 150–250 PEL while GJ 667Cc received the maximum amount of XUV radiation among them, F_XUV = 450–800 PEL.

Using the recent results by Tian (2009) that computed the exosphere properties of massive SEs, i.e., M_p ⩾ 6 M_⊕, subject to different XUV fluxes, we can estimate the exosphere radius and mass-loss rate for our three habitable SEs. Actually, since the Tian (2009) results are only available for planets with a minimum mass of M_p = 6 M_⊕, a qualitative extrapolation of the results for 10, 7, and 6 M_⊕ (see Figures 4 and 6 in his paper) shows that exobase radius and mass-loss rates are larger for less massive planets. This is particularly useful for trying to apply the Tian's results to GJ 667Cc M_p ≈ 4.5 M_⊕ and other less massive potentially habitable planets. In these cases, we will use the results by Tian (2009) to calculate a lower bound of the exobase radius and mass-loss rates.

Using the estimated XUV flux for HD 40307g (M_p ≈ 7 M_⊕) and assuming an initial CO₂-rich atmosphere, Tian's results predict that the exosphere of the planet and hence its mass-loss rate was low enough to avoid a significant early erosion of its atmosphere (see Figures 4 and 6 in his paper). This is true at least during the first 1–2 Gyr during which our magnetic protection model predicts that the planet was enshrouded by a protective magnetosphere (see Figure 9). After dynamo shutdown, the atmosphere of HD 40307g has been exposed to the direct action of the stellar wind. Assuming a stellar age of 4.5 Gyr (Tuomi et al. 2012), this effect has eroded the atmosphere for 3–4 Gyr. During this unmagnetized phase, the atmospheric mass-loss rate can be simply estimated as $\dot{M}\approx \alpha m n v_{\rm eff}$ (Zendejas et al. 2010) where α is the so-called entrainment efficiency and m, n, and v_eff are the mass, number density, and effective velocity of the stellar wind as measured at planetary distance (see Equation (2)). Using a entrainment efficiency α ∼ 0.3 (which is appropriate, for example, to describe the mass loss of the Venus atmosphere), we find that the total mass loss during the unmagnetized phase is less than 1% of a conservative estimate of the total volatile content of the planet (Tian 2009). Although our model provides only upper limits to magnetic protection and the planet could have, for example, a lighter nitrogen-rich atmosphere which is more prone to XUV-induced mass losses (Watson et al. 1981; Kulikov et al. 2006; Tian et al. 2008), this preliminary estimation suggests that HD 40307g probably still preserves a dense enough atmosphere that is able to sustain surface liquid water and hence to be actually habitable.

The case of GL 581d (M_p ≈ 6 M_⊕) is quite different. Assuming that our estimations of the XUV flux are correct, the exosphere radius predicted by Tian (2009) should be close to the actual one. In this case at times t ∼ 1 Gyr, R_exo = 1.8–2.3 R_p. However, our reference magnetosphere model predicts for this planet magnetic standoff distances R_S > 2.7 at all times. We conclude that by using our estimations GL 581d could have been protected by its intrinsic magnetic field during the critical early phases of planetary evolution and probably has preserved the critical volatiles in its atmosphere. Still, the uncertainties in the exosphere model or in the atmospheric composition and of course in the magnetic model developed here should lead to a different conclusion and further theoretical and probably observational investigation are required.

The most interesting case is that of GJ 667Cc (M_p ≈ 4.5 M_⊕). The minimum exosphere radius predicted for this planet at t ∼ 1 Gyr lies between 3.0 and 4.5 R_p while, according to our magnetosphere model, the magnetic standoff distance is R_S < 3 R_p up to 2 Gyr. Since the exosphere radius should actually be larger than that predicted with the Tian (2009) model and our magnetic model is actually optimistic, the chances that this planet was unprotected by its magnetic field in the critical first gigayear are high. But exposition does not necessarily mean a complete obliteration of the atmosphere (see, for example, the case of Venus). To evaluate the level of thermal and non-thermal obliteration of the atmosphere, we need to estimate the actual mass-loss rate. At the XUV fluxes estimated at the top of the atmosphere of this planet during the first gigayear, the minimum thermal mass-loss rate of carbon atoms from a CO₂-rich atmosphere will be larger than 2–4 × 10¹⁰ atoms cm⁻² s⁻¹ (see Figure 6 in Tian 2009). We should recall that this is actually the value for a 6 M_⊕ SE. For the actual mass of the planet, 4.5 M_⊕, the mass-loss rate could be even larger. Moreover, if, as predicted here, the exosphere is directly exposed to the stellar wind, non-thermal processes can contribute to a larger increase in the mass loss from the planetary atmosphere. At the minimum mass-loss rate, the exposed GJ 667Cc atmosphere could have lost more than ∼10⁴⁶ atoms of carbon in just ∼100 Myr and in the first gigayear the amount of carbon thermally lost to space could rise to ∼10⁴⁷ atoms. If we scale up linearly with planetary mass, the total inventory of CO₂ in the atmosphere, crust, and mantle of Venus, which is 2–3 × 10⁴⁶ molecules (see Tian 2009 and references therein), a 4.5 M_⊕ planet will have a total budget of ∼10⁴⁷ CO₂ molecules. In summary, at the minimum mass-loss rate and assuming a relatively rapid degassing of the planet, Gj 667Cc could have lost its total inventory of carbon to interplanetary space in the first couple of gigayears. Even assuming that large amounts of CO₂ are still trapped in the mantle and crust of the planet, its atmosphere should be rapidly obliterated by the stellar wind. We speculate that GJ 667Cc is a sort of "Venus-like" planet. Regardless the fact that the planet is inside the radiative HZ, it has lost its capacity to support life via a massive stellar-wind-induced loss of volatiles.

The strength of a physical model depends on the hypothesis and assumptions on which it relies. Apart from numerous albeit very common assumptions, the magnetic protection model presented here depends on three major assumptions which we discuss in the following paragraphs.

We have assumed that TPs always develop an initially metallic liquid core, irrespective of their composition and early formation history. This is not necessarily true. The formation of a metallic liquid core would depend on very complex processes and other barely known physical factors. It has been shown, for example, that under extreme water oxidation of iron the formation of a metallic core will be avoided (Elkins-Tanton & Seager 2008). In this case, silicate coreless planets will be formed. On the other hand, even if a planet is well differentiated the core could be solid from the beginning (see, e.g., VAL06). However, it should be emphasized here that our model provides only the best-case scenario of magnetic protection. Therefore, if, under the assumption of having a liquid metallic core, a planet is found to be lacking enough magnetic protection, the case when the planet is not well differentiated or when it never develops a liquid core will be even worse. In these cases, the conclusions drawn from our model will be unchanged.

The calculation of key magnetosphere properties relies on very simplistic assumptions about the complex physics behind the interaction between planetary and interplanetary magnetic fields and stellar wind. In particular, standoff distances calculated with Equation (5) assume a negligible plasma pressure inside the magnetosphere. This condition could be violated in planets with very inflated atmospheres and/or at close-in orbits. Under these extreme conditions, the magnetic standoff distance given by Equation (5) could be a poor underestimation of the actual magnetopause distance. However, the weak dependence of standoff distances and polar cap areas of the stellar wind dynamic pressure offers some idea as to the actual role that magnetospheric plasma pressure may pay in determining the size of the magnetosphere. Adding a plasma pressure term P_pl to the magnetic pressure P_mp in the left-hand side of Equation (3) is equivalent to subtracting it from the stellar wind dynamic pressure P_sw in the right-hand side. An effective stellar wind pressure $P^{\prime }_{\rm sw}=P_{\rm sw} (1-P_{\rm pl}/P_{\rm sw})$ would replace the stellar wind term in the standoff distance definition (Equation (5)). As a result a plasma pressure correcting factor (1 − P_pl/P_sw)^−1/6 will modify our estimated purely magnetic standoff distance. Even in a case where the plasma was able to exert a pressure 50% of that of the stellar wind, the standoff distance will increase by only 10%. On the other hand, in order to have a standoff distance one order of magnitude larger than that estimated using Equation (5), the plasma pressure inside the magnetosphere should amount to 99.999% of the total pressure. This is precisely what an unmagnetized planet would look like. In summary, including more realistic condition in the definition of the magnetosphere boundary will not greatly modify our results.

A final but no less important assumption in our model is that of quiet stellar wind conditions. We have only taken into account average or quiet stellar wind conditions. We have completely neglected the effects of large but transient conditions such as those produced during coronal mass ejections (CMEs). To model the effect that a steady flux/influx of CME plasma could have on close-in planets, we can modify the stellar wind pressure by maintaining the nominal velocity of the plasma but increasing the number density of wind particles by a factor of two.⁴ Taking into account that $R_S\sim P_{\rm SW}^{-1/6}$, we found that under the harsher conditions the magnetosphere radius and polar cap areas will be modified only by 10%–30% with respect to the nominal or quiet values computed here. This simplified estimation shows that our results seem to be robust in relation to uncertainties in the stellar wind pressure. However, given the complexity of the interaction between the magnetosphere and the stellar wind under active phases, a further examination of this case is required and is left to future research.

In order to study the effect that uncertainties in several critical thermal evolution parameters have in the prediction of the overall magnetic protection of potentially habitable TPs, we performed a sensitivity analysis of our model. For this purpose, we independently varied the value of six carefully chosen parameters of the model (see below) and compared the predicted dipole moment, the time of inner core formation, and the dynamo lifetime with the same values obtained using the RTEM.

We performed these comparisons for planets with five different masses: 0.7, 1.0, 3.5, 4.5, and 6.0 M_⊕ (see Figure 10). These masses correspond approximately with those of the habitable planets already discovered (see Table 2) including a hydrated Venus and present Earth. In all cases, for simplicity we assume a primordial period of rotation of P = 1 day.

Zoom In Zoom Out Reset image size

Since the dipole moment is an evolving quantity, in Figure 10 we plotted the average value of this quantity as calculated in the interval 0.7–2.0 Gyr. For times earlier than 0.7 Gyr, the stellar wind pressure is uncertain and the magnetic protection cannot be estimated (as discussed in Section 5.2, observations suggest lower stellar wind pressures at times earlier than this). For times larger than 2.0 Gyr, the flux of XUV radiation and the stellar wind pressure have decreased below the initial high levels. Although an average of the dipole moment is not phenomenologically relevant, it could be used as a proxy of the overall magnetic shielding of the planet during the harsh early phases of stellar and planetary evolution.

After studying the full set of physical parameters involved in our interior structure and thermal evolution models (see Table 1), we identified six critical parameters whose values could have noticeable effects on the results or are subject to large uncertainties. We performed an analysis of the sensitivity that the model has to the variation of the following physical parameters:

• 1.  
  The core mass fraction, CMF. This is the fraction of the planetary mass represented by the metallic core. This parameter is determined by the Fe/Si ratio of the planet that it is fixed at planetary formation or could be altered by exogenous processes (e.g., late large planetary impacts). Our reference model uses Earth's value CMF = 0.325, i.e., assumes that all planets are dominated by a silicate-rich mantle. As a comparison, Mars has a CMF = 0.23 and the value for Mercury is CMF = 0.65 (it should be noted that Mercury could have lost a significant fraction of its mantle silicates increasing the total iron fraction, probably after an early large impact). The CMF determines the size of the core and hence the thermal properties of the convective shell where the magnetic field is generated. For our sensitivity analysis, we have taken two extreme values of this parameter, CMF = 0.23 (a Mars-like core) and CMF = 0.43 (an iron-rich core). Planets with larger cores have low pressure olivine mantles and our rheological model becomes unreliable.
• 2.  
  The initial temperature contrast at the CMB, ΔT_CMB = _adbΔT_adb (see Equation (25)). This is one of the most uncertain properties in thermal evolution models. The initial core temperatures would be determined by random processes involved in the assembly and differentiation of the planet. It could vary widely from planet to planet. In order to fit the thermal history of the Earth (time of inner core formation, present size of the inner core, and magnetic field strength), we set _adb = 0.7 and applied the same value to all planetary masses. In our sensitivity analysis, we varied this parameter between two extreme values of 0.6 and 0.8. Though we are not sure that this interval is representative of planets with very different masses and formation histories, our analysis at least provides the magnitude and sign of the effect that this parameter has on the dynamo properties predicted by our thermal evolution model.
• 3.  
  High pressure viscosity rate coefficient, b (see Equation (26)). Rheological properties of silicates inside the mantle are among the most uncertain aspects of thermal evolution models. They critically determine, among other key quantities, the amount of heat that the core and mantle could transport through their respective boundary layers (see Equations (13), (20), and (22)). We found that the viscosity of the lower mantle (perovskite) is the most important source of uncertainties in our thermal evolution model. The formula used to compute viscosity at that layer (see Equation (27)) strongly depends on temperature and pressure and the parameter controlling this dependence is the "rate coefficient" b. In the RTEM, we used the value b = 12.3301 to reconstruct the thermal properties of the Earth. Using this value, all the figures reflect the strong sensitivity of the model to this parameter. To study the impact of b in the model results, we varied it in the interval 10–14.
• 4.  
  Iron thermal conductivity, k_c. This parameter controls the amount of heat coming from the core. In the RTEM, we used a value k_c = 40 W m⁻¹ K⁻¹ that fits the thermal evolution history and present magnetic field of the Earth (see Table 1). Although recent first-principles analysis suggests that values as large as 150–250 W m⁻¹ K⁻¹ could be common at Earth's core conditions (Pozzo et al. 2012), we conform here to the standard values of this parameter. Further investigations to explore values as large as those found by Pozzo et al. (2012) should be attempted. For our sensitivity analysis, we varied k_c between 35 and 70 W m⁻¹ K⁻¹, two values which are inside the typical uncertainty assumed for this property.
• 5.  
  Grüneisen parameter for iron, γ_0c. This is one of the most critical parameters of the equation of state, especially at core conditions. It strongly affects the mechanical structure, temperature profile, and phase of iron in the metallic core (for a detailed discussion on the sensitivity of interior structure models to this parameter, see, e.g., VAL06). In the RTEM, we used the reference value γ_0c = 1.36 that fits the thermal evolution history and present magnetic field of the Earth (see Table 1). Assuming different kinds of core alloys, a relatively large range of values of this parameter have been used in literature (see VAL06 and references therein). The Grüneisen parameter values have been found in the range of 1.36–2.338. Since our RTEM value is at the lower end of this range, for our sensitivity analysis we have recalculated the model for a larger value of 2.06.

Other uncertain parameters, such as the critical Rayleigh number at the CMB, Ra_c, which is also varied to study the sensitivity of thermal evolution models (E. Gaidos 2011, private communication), were also studied. We did not find significant sensitivity of our model to the variation of those parameters.

In Figure 10 we depict the relative variation in the aforementioned magnetosphere and dynamo properties when each of the previously described parameters were varied independently.

We found that planetary composition (CMF) and mantle viscosity are responsible for the largest uncertainties in the predicted magnetic properties of the planet. Planets with small metallic cores have on average low magnetic dipole moments (squares in the first column of the upper panel). This is mainly due to a geometrical effect. The total heat produced by the core and hence the magnetic field strength at the core surface is of the same order as Fe-poor and Fe-rich planets. However, a small core also means a lower magnetic dipole moment, i.e., ${\cal M}\sim R_c^3$. Planets with lower iron contents also have small and hot cores and therefore the solid inner core formation and shutting down the dynamo are slightly delayed (squares in the middle and lowest panel of Figure 10).

Viscosity dependence on pressure and temperature, as quantified by the parameter b, has the opposite effect on planetary magnetic properties than CMF at least for EPs. A low viscosity, lower mantle will favor the extraction of heat from the metallic core increasing the available convective energy for dynamo action. On the other hand, a viscous lower mantle will delay the formation of a solid inner core and extend the lifetime of the dynamo (middle and lowest panel in Figure 10).

The effect of the Grüneisen parameter at core conditions is negligible, at least in respect to the magnetic field strength and lifetime (upper and lower panels) which are the most critical properties affecting planetary magnetic protection. Only the time of the inner core formation is strongly affected by changes in this parameter. In planets smaller than Earth, inner core solidification can be delayed up to three times the reference value. On the other hand, with a larger Grüneisen parameter, the inner cores of EPs could become solid very early in their thermal histories, even almost from the beginning. This is a result of the interplay between the resulting evolution of the thermal profile and the solidus.

The effect of the initial temperature contrast across the CMB, quantified by the parameter _adb, goes in the same direction as viscosity. The reasons for this behavior are, however, far more complex. A larger initial temperature contrast across the CMB also implies a larger initial temperature at the core center. Although a hotter core also produces a larger amount of available convective energy, the time required for iron to reach the solidification temperature is also larger. The dynamo of planets with M_p < 2 M_⊕ and hot cores (large _adb) is weaker than that of more massive planets during the critical first couple of gigayears where the average is calculated. Planets with colder cores develop solid inner cores almost from the beginning and the release of latent and gravitational energy feeds stronger dynamos.

The case of more massive planets, M_p > 2 M_⊕, where the condition for an inner core formation is never reached during the dynamo lifetime, is different. In this case, planets with hot cores (large lower mantle viscosities or high temperature contrasts along the CMB) produce large amounts of available convective energy. A larger convective power will produce a larger value of the local Rossby number. Thus, massive planets with hot cores also have multipolar dynamos and hence lower dipole magnetic moments and a reduced magnetic protection.

Thermal conductivity affects the results of the thermal evolution model less. Besides the case of massive planets where differences on the order of 10%–30% in the magnetic properties are observed when k_c varies between its extremes, the magnetic properties calculated with our RTEM seem very robust in comparison with variations of these two properties. However, it should be mentioned that this result applies only when a standard value of k_c is assumed. Further investigations to explore the recent findings (Pozzo et al. 2012) concerning the possibility that k_c could be larger by factors as large as 2–3 should be attempted.

In summary, despite the existence of a natural sensitivity of our simplified thermal evolution model to uncertainties on their free parameters, the results presented in this paper seem to be correct at least on the order of magnitude. Moreover, since the standoff distance and polar cap areas, which are the actual proxies to magnetic protection, goes as ${\cal M}^{1/3}$, a one order of magnitude estimation of ${\cal M}$ will give us an estimation of the level of magnetic protection by a factor of around two.

To clarify this point, let us consider the case of GL 581d. If we assume, for example, that its iron content is much less than that of the Earth (as was assumed in the RTEM model), but the rest of critical thermal properties are essentially the same, the average standoff distance (polar cap area) at the critical first gigayear will be off by only ∼30% with respect to the prediction depicted in Figure 9. More interesting is to note that the uncertainties due to the unknown period of rotation (shaded area in Figure 9) seem to be much larger than those coming from the uncertainties in the thermal evolution model.

Validating or improving thermal evolution models and dynamo scaling laws for the case of SEs presently represents a huge observational challenge. The available sensitivity of our best instruments on Earth and in space is rather insufficient. New and/or improved instruments and observational techniques will be required to detect, catalog, and compare the thermal and magnetic properties of low-mass planets in the future. However, the importance that the detection and characterization of the magnetic properties of potentially habitable planets discovered in the future in order to assess their true habitability clearly justifies the effort.

The first goal seems to be the direct or indirect detection of SE magnetospheres. Four methods, some of which are already used in our own solar system, could be devised to achieve this goal: (1) the detection of radio waves coming from synchrotron and cyclotron radiation produced by plasma trapped in the magnetosphere, (2) the detection of a bow shock or a tail of ions produced by the interaction of the planetary atmosphere and magnetosphere with the stellar wind or the interplanetary plasma, (3) the detection of planetary auroras, and (4) spectroscopic observations of a non-equilibrium atmospheric chemistry induced by a high flux of CR (this is actually a negative detection of a magnetosphere).

The first (radio emission) and second methods (bow shock or tail) have already been studied in detail (Bastian et al. 2000; Farrell et al. 2004; Grießmeier et al. 2007b; Lazio et al. 2009; Vidotto et al. 2010, 2011). Its reliability, at least for the case of planets with intense magnetic fields or placed very close to their host star, has been already tested.

If synchrotron or cyclotron radiation coming from the magnetosphere of TPs could be detected, the power and spectra of the radiation could be used to measure the magnetic field strength. However, even with the most sensitive instruments, e.g., the Low Frequency Array or the Long Wavelength Array, the expected power and spectra are several orders of magnitude below the threshold of detection. Powers as large as 10³–10⁵ times the Jupiter radio emission and frequencies in the range of tens of MHz are required for the present detection of synchrotron radio emission in planetary magnetospheres (see, e.g., Grießmeier et al. 2007b). The magnetic field intensities and expected frequencies produced in SEs magnetospheres are several orders of magnitude lower than these thresholds and are not likely to be detected in the near future.

It has been shown recently that measurements of the asymmetry in the ingress and egress of transiting planets can be used to detect the presence of a bow shock or a tail of plasma around the planet. Vidotto et al. (2010, 2011) used this phenomenon to constrain the magnetic properties of Wasp-12b. The formation of a detectable bow shock depends, among other factors, on the relative velocity between the planet and the shocked plasma. Close-in planets with strong enough magnetic fields (this is precisely the case of Wasp-12b) can easily develop UV-opaque bow shocks and allow reliable detection. However, low-mass planets with relatively weak magnetic fields, such as those predicted with our models, hardly produce a detectable bow shock. It has been estimated that magnetopause fields in the range of several gauss (among other important factors such as large enough projected magnetospheric size) should be required to have a detectable signal of a bow shock (A. A. Vidotto 2012, private communication). Our habitable SEs have magnetopause fields of the order of a few microgauss (see Figure 7). The case of an ion tail coming from a weakly magnetized planet has received less attention and probably could offer better chances for a future indirect detection of the magnetic environment of low-mass planets.

Finally, the detection of far-UV (FUV) or X-ray emission from planetary auroras can also be used as a tool to directly and indirectly study planetary magnetospheres. Planetary auroras with intensities as high as 10²–10³ times larger than that of the Earth are expected in close-in giant planets subject to the effects of CMEs from its host star (Cohen et al. 2011). If we estimate that a typical Earth aurora has an intensity of 1 kR (Neudegg et al. 2001; being 1 R ∼10⁻¹¹ photons m⁻² s⁻¹ sr⁻¹) and assuming that 10% of a close-in Jupiter-like planet is covered by auroras producing FUV photons around 130 nm, the total emitted power from these planets will be ∼10¹³ W. If we assume that this is the present threshold for exoplanetary aurora detection, even under strong stellar wind conditions and distances typical of the HZ, the total FUV power produced by auroras in the polar cap of EPs could be only 10⁵ W, which is eight orders of magnitude less than the present detection threshold. In addition, the FUV radiation should be detected against an intense UV background likely coming from a young and active low-mass star. If we can find ways to overpass these difficulties, the observation of the FUV and X-ray emission and its variability from auroras in potentially habitable SEs could be used as powerful probes of the magnetic environment around the planet.