{'fraction_non_alphanumeric': 0.027901319393506097, 'fraction_numerical': 0.017743567840595723, 'mean_word_length': 4.5260133401744485, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 4, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Influence of Thermal Evolution in the Magnetic Protection of Terrestrial Planets

Magnetic protection of potentially habitable planets plays a central role in determining their actual habitability and/or the chances of detecting atmospheric biosignatures. We develop here a thermal evolution model of potentially habitable Earth-like planets and super-Earths. Using up-to-date dynamo scaling laws we predict the properties of core dynamo magnetic fields and study the influence of thermal evolution on their properties. The level of magnetic protection of tidally locked and unlocked planets is estimated by combining simplified models of the planetary magnetosphere and a phenomenological description of the stellar wind. Thermal evolution introduces a strong dependence of magnetic protection on planetary mass and rotation rate. Tidally locked terrestrial planets with an Earth-like composition would have early dayside magnetospause distances between 1.5 and 4.0 Rp, larger than previously estimated. Unlocked planets with periods of rotation ~1 day are protected by magnetospheres extending between 3 and 8 Rp. Our results are robust against variations in planetary bulk composition and uncertainties in other critical model parameters. For illustration purposes the thermal evolution and magnetic protection of the potentially habitable super-Earths GL 581d, GJ 667Cc and HD 40307g were also studied. Assuming an Earth-like composition we found that the dynamos of these planets are already extinct or close to being shut down. While GL 581d is the best protected, the protection of HD 40307g cannot be reliably estimated. GJ 667Cc, even under optimistic conditions, seems to be severely exposed to the stellar wind and, under the conditions of our model, has probably suffered massive atmospheric losses.


INTRODUCTION
The discovery of extrasolar habitable planets is one of the most ambitious challenges in exoplanetary research. At the time of writing, there are almost 861 confirmed exoplanets 5 including 61 classified as Earth-like planets (EPs, M ∼ 1 M ⊕ ) and super-Earths (SEs, M ∼ 1 − 10 M ⊕ , Valencia et al. 2006 hereafter VAL06). Among these low mass planets there are three confirmed SEs, GJ 667Cc (Bonfils et al. 2011), GL 581d (Udry et al. 2007;Mayor et al. 2009), and HD 40307g (Tuomi et al. 2012), and tens of Kepler candidates (Borucki et al. 2011;Batalha et al. 2012) that are close or inside the habitable zone (HZ) of their host stars (see e.g. Selsis et al. 2007;Pepe et al. 2011;Kaltenegger et al. 2011). If we include the possibility that giant exoplanets could harbour habitable exomoons, the number of the already discovered potentially habitable planetary environments beyond the Solar System could be rised to several tens (Underwood et al. 2003;Kaltenegger 2010). Moreover, the existence of a plethora of other terrestrial planets 1 jzuluaga@fisica.udea.edu.co 2 sbustama@pegasus.udea.edu.co 3 p.cuartas@fisica.udea.edu.co 4 jhhoyos@udem.edu.co 5 For updates, please refer to http://exoplanet.eu (TPs) and exomoons in the Galaxy is rapidly gaining evidence (Borucki et al. 2011;Catanzarite & Shao 2011;Bonfils et al. 2011;Kipping et al. 2012), and the chances that a large number of potentially habitable extrasolar bodies could be discovered in the near future are very large.
The question of which properties a planetary environment needs in order to allow the appearance, evolution and diversification of life has been extensively studied (for recent reviews see Lammer et al. 2009 andKasting 2010). Two basic and complementary physical conditions must be fulfilled: the presence of an atmosphere and the existence of liquid water on the surface (Kasting et al. 1993). However, the fulfilment of these basic conditions depends on many complex and diverse endogenous and exogenous factors (for a comprehensive enumeration of these factors see e.g. Ward & Brownlee 2000 The existence and long-term stability of an intense planetary magnetic field (PMF) is one of these relevant factors (see e.g. Grießmeier et al. 2010 and references therein). It has been shown that a strong enough PMF could protect the atmosphere of potentially habitable planets, especially its valuable content of water and other volatiles, against the erosive action of the stellar wind (Lammer et al. 2003Khodachenko et al. 2007;Chaufray et al. 2007). Planetary magnetospheres also act as shields against the potentially harmful effects that the stellar high energy particles and galactic cosmic rays (CR) produce in the life-forms evolving on the planetary surface (see e.g. Grießmeier et al. 2005). Even in the case that life could arise and evolve on unmagnetized planets, the detection of atmospheric biosignatures would be also affected by a higher flux of high energy particles including CR, especially if the planet is close to very active M-dwarfs (dM) (Grenfell et al. 2007;Segura et al. 2010).
It has been recently predicted that most of the TPs in our Galaxy could be found around dM stars (Boss 2006;Mayor & Udry 2008;Scalo et al. 2007;Rauer et al. 2011;Bonfils et al. 2011). Actually ∼ 20% of the presently confirmed super-Earths belong to planetary systems around stars of this type. Planets in the HZ of low mass stars (M ⋆ 0.6M ⊙ ) would be tidally locked (Joshi et al. 1997;Heller et al. 2011), a condition that poses serious limitations to their potential habitability (see e.g. Kite et al. 2011 and references therein). Tidally locked planets inside the HZ of dMs have periods in the range of 5 − 100 days, a condition that has commonly been associated with the almost complete lack of a protective magnetic field ). However, the relation between rotation and PMF properties, that is critical at assessing the magnetic protection of slowly rotating planets, is more complex than previously thought. In particular, a detailed knowledge of the thermal evolution of the planet is required to predict not only the intensity but also the regime (dipolar or multipolar) of the PMF for a given planetary mass and rotation rate (Zuluaga & Cuartas 2012).
Although several authors have extensively studied the protection that intrinsic PMF would provide to extrasolar planets (Grießmeier et al. 2005;Khodachenko et al. 2007;Lammer et al. 2007;Grießmeier et al. 2009Grießmeier et al. , 2010, all have disregarded the influence that thermal evolution has in the evolution of planetary magnetic properties. They have used also outdated dynamo scaling-laws that have been recently revised (see Christensen 2010 and references therein). The role of rotation in determining the PMF properties that is critical in assessing the case of tidally locked planets has also been overlooked (Zuluaga & Cuartas 2012).
We develop here a comprenhensive model for the evolution of the magnetic protection of potentially habitable TPs around GKM main sequence stars. To achieve this goal we integrate in a single framework a parametrized thermal evolution model based on the most recent advances in the field (Gaidos et al. 2010;Tachinami et al. 2011;Stamenković et al. 2012), up-to-date dynamo scaling-laws (Christensen 2010;Zuluaga & Cuartas 2012) and phenomenological models for the evolution of the stellar wind and planet-star magnetic interaction (Grießmeier et al. 2010). Our model is aimed at: 1) understanding the influence of thermal evolution in the magnetic protection of TPs, 2) assesing the role of low rotation periods in the evolution of the magnetic protection of tidally-locked habitable planets, 3) placing more realistic constraints on the magnetic properties of potentially habitable TPs suitable for future studies of atmospheric mass-loss or the CR effect on the atmospheric chemistry or on life itself and 4) estimat-ing by the first time the magnetic properties of already discovered super-Earths in the HZ of their host stars.
Our work is a step-forward in the understanding of planetary magnetic protection because it puts together in a single model the evolution of the magnetic properties of the planet and hence its dependence on planetary mass and composition and the role of the planet-star interaction into determining the resulting level of magnetic protection. Previous models of the former (thermal evolution and intrinsic magnetic properties) did not consider the interaction of the planetary magnetic field with the evolving stellar wind which is finally the factor that determine the actual level of planetary magnetic protection. On the other hand previous attempts to study the planetstar interaction overlooked the non-trivial dependence of intrinsic magnetic properties on planetary thermal evolution and hence on planetary mass and rotation rate. Additionally, and for the first time, we are attempting here to calculate the magnetic properties of the already discovered potentially habitable super-Earths, GJ 667Cc, Gl 581d and HD 40307g.
This paper is organized as follows. Section 2 is aimed at introducing the properties of planetary magnetospheres we should estimate in order to evaluate the level of magnetic protection of a potentially habitable terrestrial planet. Once those properties are expressed in terms of two basic physical quantities, the planetary magnetic dipole moment and the pressure of the stellar wind, we proceed to describe how those quantities can be estimated by modelling the thermal evolution of the planet (section 3), scaling the dynamo properties from the planetary thermal and rotational properties (section 4), and modelling the interaction between the star and the planet (section 5). In section 6 we apply our model to evaluate the level of magnetic protection of hypothetical potentially habitable TPs as well as the already discovered habitable super-Earths. We also present there the results of a numerical analysis aimed at evaluating the sensitivity of our model to uncertainties in the composition of the planet and to other critical parameters of the model. In section 7 we discuss the limitations of our model, present an example of the way in which our results could be applied to estimate the mass-loss rate from already and yet to be discovered potentially habitable TPs and discuss the observational prospects to validate or improve the model. Finally a summary and several conclusions drawn from this research are presented in section 8.

MAGNETOSPHERE
The interaction between the PMF, the interplanetary magnetic field (IMF) and the stellar wind creates a magnetic cavity around the planet known as the magnetosphere. Although magnetospheres are very complex systems, their global properties are continuous functions of only two physical variables (Siscoe & Christopher 1975): the magnetic dipole moment of the planet, M, and the dynamic pressure of the stellar wind, P sw . Dipole moment is defined in the multipolar expansion of the magnetic field strength: where B p (r) is the angular-averaged PMF strength measured at a distance r from the planet center and µ 0 = 4π × 10 −7 H/m is the vacuum permeability. In the rest of the paper we will drop-off the higher order terms in 1/r (multipolar terms) and focus on the dipolar component of the field B dip p which is explicitly given by the first term of the right side in eq. (1).
The dynamic pressure of the stellar wind is given by where m and n are the typical mass of a wind particle (mostly protons) and its number density, respectively.
Here v eff = (v 2 sw + v 2 p ) 1/2 is the effective velocity of the stellar wind as measured in the reference frame of the planet whose orbital velocity is v p . T is the local temperature of the wind plasma and k B = 1.38 × 10 −23 j/K is the Boltzmann constant.
There are three basic properties of planetary magnetospheres we are interested in: 1) the maximum magnetopause field intensity B mp , which is a proxy of the flux of high energy particles entering into the magnetospheric cavity, 2) the standoff or stagnation radius, R S , a measure of the size of the dayside magnetosphere, and 3) the area of the polar cap A pc that measures the total area of the planetary atmosphere exposed to open field lines through which particles can escape to interplanetary space. The value of these quantities provides information about the level of exposure that a habitable planet has to the erosive effects of stellar wind and the potentially harmful effects of the CR.
The maximum value of the magnetopause field intensity B mp is estimated from the balance between the magnetic pressure P mp = B 2 mp /(2µ 0 ) and the dynamic stellar wind pressure P sw (eq. (2)), Here we are assuming that the pressure exerted by the plasma inside the magnetospheric cavity is negligible (see discussion below).
Although magnetopause fields arise from very complex processes (Chapman-Ferraro and other complex currents at the magnetosphere boundary), in simplified models B mp is assumed proportional to the PMF intensity B p as measured at the substellar point r = R S (Mead 1964;Voigt 1995), f 0 is a numerical enhancement factor of order 1 that can be estimated numerically. We are assuming here that the dipolar component of the intrinsic field (first term in the r.h.s. of eq. (1)) dominates at magnetopause distances even in slightly dipolar PMF.
Combining equations (3) and (4) we estimate the standoff distance: sw that can be expressed in terms of the present dipole moment of the Earth M ⊕ = 7.768 × 10 22 A m 2 and the average dynamic pressure of the solar wind as measured at the orbit of our planet P sw⊙ = 2.24 × 10 −9 (Stacey 1992;Grießmeier et al. 2005): It is important to stress that the value of R S estimated with eq. (5) assumes a negligible value of the plasma pressure inside the magentospheric cavity. This approximation is valid if at least one of these conditions is fulfilled: 1) the planetary magnetic field is very intense, 2) the dynamic pressure of the stellar wind is small, or 3) the planetary atmosphere is not too bloated by the XUV radiation. In the case when any of these conditions are fulfilled we will refer to R S as given by eq. (5) as the magnetic standoff distance which is an underestimation of the actual size of the magnetosphere.
Last but not least we are interested in evaluating the area of the polar cap. This is the region in the magnetosphere where magnetic field lines could be open into the interplanetary space or to the magnetotail region. Siscoe & Chen (1975) have shown that the area of the polar cap A pc scales with dipole moment and dynamic pressure as: Here we have normalized the polar cap area with the total area of the atmosphere 4πR 2 p , and assumed the atmosphere has a scale-height much smaller than planetary radius R p . In order to model the evolution of these three key magnetosphere properties we need to estimate the surface dipolar component of the PMF B dip p (R p ) (from which we can obtain the dipole moment M), the average number density n, velocity v eff and temperature T of the stellar wind (which are required to predict the dynamic pressure P sw ). These quantities depend in general on time and also on different planetary and stellar properties. In the following sections we describe our model for the calculation of the evolving values of these fundamental quantities.

THERMAL AND DYNAMO EVOLUTION
We assume here that the main source of a global PMF in TPs is the action of a dynamo powered by convection in a liquid metallic core (Stevenson 1983(Stevenson , 2003. This assumption is reasonable since the Moon and all rocky planets in the solar system, regardless of their different origins and compositions, seem to have in the present or to have had in the past an iron core dynamo (see e.g. Stevenson 2010). Other potential sources of PMFs such as body currents induced by the stellar magnetic field or dynamo action in a mantle of ice, water or magma are not considered here but left for future research.
The properties and evolution of a core dynamo will depend on the internal structure and thermal history of the planet. Thermal evolution of TPs, specially the Earth itself, has been studied for decades (for a recent review see Nimmo 2009). A diversity of thermal evolution mod-els for planets larger than the Earth have recently appeared in literature (Papuc & Davies 2008;Gaidos et al. 2010;Tachinami et al. 2011;Driscoll & Olson 2011;Stamenković et al. 2011). But the lack of observational evidence against which we can compare the predictions of these models has left too much room for uncertainties especially regarding mantle rheology, core composition and thermodynamic properties. Albeit these fundamental limitations, a global picture of the thermal history of super-Earths has started to arise. Here we follow the lines of Labrosse 2003 andGaidos et al. 2010 developing a parametrized thermal evolution model which combines a simplified model of the interior structure and a parametrized description of the core and mantle rheology.
Our model includes several distinctive characteristics in comparison to previous ones. The most important one is an up-to-date treatment of mantle rheology. For that purpose we use two different formulae to compute the viscosity of the upper and lower mantle. By lower mantle we understand here the region of the mantle close to the core mantle boundary (CMB). This is the hottest part of the mantle. The upper mantle is the outer cold part of this layer. It is customary to describe both regions with the same rheology albeit their very different mineralogical compositions. Additionally thermal and density profiles in the mantle following the same prescription as in the core. We also use a different ansatz to assign initial values to lower mantle temperature and to the temperature contrast across the CMB, two of the most uncertain quantities in thermal evolution models. Using our ansatz we avoid assigning arbitrary initial values to these critical parameters but more importantly we are able to find a unified method to set the value of these temperatures in planets with very different masses. It is also important to notice that in other models these temperatures were set by hand or were treated as free parameters in the model Four key properties should be predicted by any thermal evolution model in order to calculate the magnetic properties of a planet: 1) the total available convective power Q conv , providing the energy required for magnetic field amplification through dynamo action; 2) the radius of the solid inner core R ic and from there the height D ≈ R c − R ic of the convecting shell where the dynamo action takes place (R c is the radius of the core); 3) the time of inner core formation t ic and 4) the total dynamo life-time t dyn .
In order to calculate these quantities we solve simply parametrized energy and entropy equations of balance describing the flux of heat and entropy in the planetary core and mantle.
As stated before our model is based on the interior structure model by VAL06) and in thermal evolution models previously developed by Schubert et al. 1979, Stevenson 1983, Nimmo & Stevenson 2000, Labrosse et al. 2001, Labrosse 2003, Gubbins et al. 2003, Gubbins et al. 2004, Aubert et al. 2009,Gaidos et al. 2010,Stamenković et al. 2011. For a detailed description of the fundamental physics behind the thermal evolution model developed here please refer to these earlier studies.

Interior structure
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 solving simultaneously the continuity, Adams-Williamson and hydrostatic equillibrium equations (eqs. (1)-(4) in VAL06). For all planetary masses we assume boundary conditions, ρ(r = R p ) = 4000 kg m −3 and P (r = R p ) = 0 Pa. For each planetary mass and core mass fraction CMF= M core /M p we use a RK4 integrator and a shooting method to compute consistently 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 of that we assume a mantle completely made of perovskite-postperovskite (ppv). This is the reason why we take ρ(r = R p ) = 4000 kg m −3 instead of the more realistic value of 3000 kg m −3 . 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 super-Earths , i.e. 100 − 1000 GPa for M p = 1 − 10M ⊕ (VAL06 and 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 constant during the whole thermal evolution of the planet. We have verified that this assumption does not change significantly the mechanical properties inside the planet at least not at a level affecting the thermal evolution itself.
Inside the core we update continuously 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 at 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.

3.2.
Core thermal evolution 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): 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 are negligible at the typical conditions inside super-Earths (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 turn-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: 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 : 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): where D c = 3c p /2παρ 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 eq. (9): 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): 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), 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 reference therein). This temperature could be computed in general as a weighted average of the temperatures at the boundaries, 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 parametrized by VAL06: Here δ(ρ) ≈ 1/3 and γ is an effective Grüneisen parameter that is assumed for simplicity constant. To integrate this equation we use the numerical density profile provided by the interior model and the reference values ρ 0 = 8300 kg/m 3 (pure iron) and τ 0 = 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 eqs. (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), Here ∆ is the ratio between the gradient of the solidus (eq. (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 ⋆ (eq. (27) in Gaidos et al. 2010) and the temperature at the core surface is increase to T ⋆ (eq. (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, where the total entropy Φ is computed from the Entropy balance (eq. (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 eq. (8), Q conv gets also negative and convection is not 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 life-time t dyn .

Mantle thermal evolution
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: 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 (eq. (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 well the Earth properties. 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 analogous expression to eq. (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 right below the surface boundary layer (see figure 1).
Assuming an adiabatical temperature profile in the mantle we can also write: in analogy to the core temperature profile, eq. (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 2 m = L 2 m /γ Labrosse (2003). In our simplified model we take the values of the density at the boundaries of the mantle and obtain analytically 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 the presence or not of mobile lids play an important role into 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 , where ∆T m = T m − T s is the temperature contrast across the SB and T s the surface temperature. Since we are studying the thermal evolution of habitable planets we assume in all cases T s = 290 K. Planetary interior structure and thermal evolution are not too sensitive to surface temperature. We have verified that results are nearly the same for 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), 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: Here Γ ≡ −∂ ln η m /∂T m measures the viscosity dependence on temperature evaluated at the average mantle pressure.
With all this elements at hand the energy balance at eq. (19) is finally transformed into an ordinary differential equation for the upper mantle temperature T m (t), where C m is the bulk heat capacity of the upper mantle which is calculated with an analogous expression to eq. (12)

Initial conditions
In order to solve the coupled differential eqs. (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 , Here θ = 1700 K is a potential temperature which is assumed the same for all planetary masses (Stamenković et al. 2011). Using T m (t = 0) and the adiabatic temperature profile we can get the initial lower mantle temperature T l (t = 0).
The initial value of core temperature T c (t = 0) is one of the most uncertain parameters in thermal evolution models. Although nobody knows its actual value or its dependence on the formation history and planetary mass, 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 3 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 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: 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 criterium 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 .

Rheological model
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 sillicates and iron at high pressures and temperatures. A detailed discussion on this important topic is out of the scope of this paper. An up-to-date discussion and analysis of the dependence on pressure and temperature of viscosity in super-Earths and its influence in thermal evolution can be found in the recent works by Tachinami et al. 2011 andStamenković 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), the grain size and m the growing exponent, A and b are free parameters, D 0 is the pre-exponential diffusion coefficient and 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: All the parameters used in the visocisty 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 eq. (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 mantle avoids an 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 Arrheniustype 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 highpressures 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 as 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 parametrization given by Tachinami et al. (2011): (27) whereǫ 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 −6 m 3 mol −1 . Using the formula in eq. (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 are 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 − 6M ⊕ .

PLANETARY MAGNETIC FIELD
In recent years improved numerical experiments have constrained the full set of possible scaling laws used to predict the properties of planetary and stellar convectiondriven dynamos (see Christensen 2010 and references therein). It has been found that in a wide range of physical conditions the global properties of a planetary dynamo can be expressed in terms of simple power-law functions of the total convective power Q conv and the size of the convective region.
One of the most important results of power-based scaling laws is the fact that the volume averaged magnetic field intensity B 2 rms = (1/V ) B 2 dV does not depend on the rotation rate of the planet (eq. (6) in Zuluaga & Cuartas 2012), Here C Brms is a fitting constant obtained from numerical dynamo experiments and its value is different in the case of dipolar dominated dynamos, C dip Brms = 0.24, and multipolar dynamos, C mul are the average density, height and volume of the convective shell.
The dipolar field intensity at the planetary surface, and hence the dipole moment of the PMF, can be estimated if we have information about the power spectrum of the magnetic field at the core surface. Although we cannot predict the relative contribution of each mode to the total core field strength, numerical dynamos exhibit an interesting property: there is a scalable dimensionless quantity, the local Rossby number Ro * l , that could be used to distinguish dipolar dominated from multipolar dynamos. The scaling relation for Ro * l is (eq. (5) in Zuluaga & Cuartas 2012): (29) Here C Rol = 0.67 is a fitting constant and P is the period of rotation. It has been found that dipolar dominated fields arise systematically when dynamos have Ro * l < 0.1. Multipolar fields arise in dynamos with values of the local Rossby number close to and larger than this critical value. From eq. (29) we see that in general fast rotating dynamos (low P ) have dipolar dominated core fields while slowly rotating ones (large P ) produce multipolar fields and hence fields with a much lower dipole moment.
It is important to stress that the almost independence of B rms on rotation rate, together with the role that rotation has in the determination of the core field regime, implies that even very slowly rotating planets could have a magnetic energy budget of comparable sized than rapidly rotating planets with similar size and thermal histories. In the former case the magnetic energy will be redistributed among other multipolar modes rendering the core field more complex in space and probably also in time. Together all these facts introduce a non-trivial dependence of dipole moment on rotation rate very different from that obtained with the traditional scaling laws used in previous works (see e.g. Grießmeier et al. 2004 andKhodachenko et al. 2007). Here we want to emphasize a property that was also previously overlooked. Multipolar dominated dynamos produces magnetic fields that decays more rapidly with distance than dipolar fields and so it is expected that a planet with a multipolar magnetic field will be less protected than those having strongly dipolar fields.
Using the value of B rms and Ro * l we can compute the maximum dipolar component of the field at core surface. For this purpose we use an upper bound to the dipolarity fraction f dip (the ratio of the dipolar component to the total field strength at core surface). Dipolar dominated dynamos have by definition f dip ≤ f max dip = 1.0. The case of reversing dipolar and multipolar dynamos is more complex. Numerical dynamo experiments show that multipolar dynamos have Ro * l 0.1 and f dip f max dip = 0.35. However to avoid inhomogeneities in the transition region around Ro * l ≈ 0.1 we calculate a maximum dipolarity fraction through a "soft step function", To connect this ratio to the volumetric averaged magnetic field B rms we use the volumetric dipolarity fraction b dip that it is found, as shown by numerical experiments, conveniently related with the maxium value of where c bdip ≈ 2.5 is again a fitting constant. It is important to notice here that the exponent 11/10 is the ratio of the smallest integers close to the numerical value of the fitting exponent (see figure 1 in Zuluaga & Cuartas 2012 Radiactive heat correction 1.253 -  Finally by combining eqs. (28)-(30) we can compute an upper bound to the dipolar component of the field at the CMB: The surface dipolar field strength is estimated using, and finally the total dipole moment is calculated using eq. (1) for r = R p .
It should be emphasized that the surface magnetic field intensity determined using eq. (32) overestimates the PMF dipolar component. The actual field could be much more complex spatially and the dipolar component could be lower. As a consequence our model can only predict the maximum level of protection that a given planet could have from a dynamo-generated intrinsic PMF.
The results of applying the RTEM to calculate the properties of the magnetic field of TPs in the mass range 0.5-4.0 M ⊕ using the scaling laws in eqs. (28), (29) and (31) are summarized in figures 3 and 4. In figure 3 we show the local Rossby number, the maximum dipolar field intensity and the dipole moment as a function of time computed for planets with different mass and two different periods of rotation (P = 1 day and P = 2 days). This figure shows the effect that rotation has on the evolution of dynamo geometry and hence in the maximum attainable dipolar field intensity at the planetary surface. In figure 4 we have summarized in mass-period (M-P) diagrams (Zuluaga & Cuartas 2012) the evolution of the dipole moment for planets with long-lived dynamos. We see that for periods lower than 1 day and larger than 5-7 days the dipole moment is nearly independent of rotation. Slowly rotating planets have a non-negligible dipole moment which is systematically larger for more massive planets.

PLANET-STAR INTERACTION
The PMF properties constrained using the thermal evolution model and the dynamo-scaling laws are not enough to evaluate the level of magnetic protection of a potentially habitable TP. We also need to estimate also the magnetosphere and stellar properties (stellar wind and luminosity) as a function of time in order to assess properly the level of star-planet interaction.
Since the model developed in previous sections provides only the maximum intensity of the PMF, we will be interested here into constraint the magnetopshere and stellar properties from below, i.e. to find the lower level of "stellar aggression" for a given star-planet configura-tion. Combining upper bounds of PMF properties and lower bounds for the star-planet interaction will produce an overestimation of the overall magnetic protection of a planet. If under this model a given star-planet configuration is not suitable to provide enough magnetic protection to the planet, the actual case should be much worse. But if, on the other hand, our upper-limit approach predicts a high level of magnetic protection, the actual case could still be that of an unprotected planet. Therefore our model is capable at predicting which planets will be unprotected but less able when predicting which ones will be actually protected.

The Habitable Zone (HZ) and tidally locking limits
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 starplanet distance (distance to the HZ) and 3) the conmensurability of planetary rotation and orbital period (tidal locking). These properties should be properly modelled in order to assess the degree of star-planet interaction which are critical at 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 in assessing the actual properties of newly discovered exoplanets. The case of low mass main sequence stars (GKM) are particularly important in 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 by 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 to 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 eq. (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.
The orbital and rotational properties of planets at close-in orbits are strongly affected by the gravitational and tidal interaction 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 , 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 estimating in general 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): 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 hours (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 already discovered potentially habitable super-Earths, GL 581d, GJ 667Cc and HD 40307g, are also highlighted in this figure.

Stellar wind
The stellar wind and cosmic rays pose the highest risks for a magnetically unprotected potentially habitable terrestrial planet. The dynamic pressure of the wind is able to obliterate an exposed atmosphere, especially druing the early phase of stellar evolution (Lammer et al. 2003), and energetic stellar cosmic rays 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 developed originally by Parker 1958 that describes the wind as a non-magnetized, isothermal and axially symmetric flux of particles (hereafter the Parker's model) and the more detailed albeit simpler magneto-hydrodynamic model developed originally 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 describes reliably the properties of the stellar winds in the case of stars with periods of rotation of the same order of 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 almost by a factor of 2 (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 the Parker's wind equation (Parker 1958): where u = v/v c and ρ = d/d c are the velocity and distance normalized with respect to v c = k B T /m and d c = GM ⋆ m/(4k B T ) 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 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: HereṀ ⋆ 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Ṁ ⋆ .
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 andLammer 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 Here α v = −0.43, α n = −1.86 ± 0.6 and τ = 25.6 Myr (Grießmeier et al. 2009). The parameters v 0 = 3971 km/s and n 0 = 1.04 × 10 11 m −3 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 6 m −3 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Ṁ ⋆ (t) in the 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 here this procedure. 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 1AU , is calculated using equation (37). Replacing this velocity in the Parker's wind equation for d = 1 AU, we find numerically the temperature of the Corona T (t). This parameter is enough to provide us 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 eqs. (37) and (38) the mass-loss rate for the SunṀ ⊙ at time t and d = 1 AU can be obtained: Assuming that the mass-loss rate scales-up simply with the stellar surface area, i.e.Ṁ ⋆ (t) =Ṁ ⊙ (t)(R ⋆ /R ⊙ ) 2 , the value ofṀ ⋆ can finally be estimated.
Using v(t, d, M ⋆ ) andṀ ⋆ 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 ⋆ ) 2 inside the HZ of four different stars as computed using the procedure described before is plotted in figure 6.
It is important to stress here that for stellar ages t 0.7 Gyr the semiempirical formulae in eqs. (37) and (38) are not longer reliable (Grießmeier et al. 2007a). These equations are based in the empirical relationship observed between the X-ray surface flux and the massloss rateṀ ⋆ (Wood et al. 2002(Wood et al. , 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, 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 are 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.

RESULTS
Using the results of our RTEM, the power-based scaling laws for dynamo properties, and the properties of the stellar insolation and stellar wind pressure, we have calculated the magnetosphere properties of earth-like planets and super-Earths in the HZ of different mainsequence stars. We have performed these calculations for hypothetical TPs in the mass-range 0.5-6 M ⊕ and for the already discovered potentially habitable planets GL 581d, GJ 667Cc and HD 40307g (see table 2). The case of the Earth and an habitable Venus has also been studied for references purposes.
To include the effect of rotation in the properties of the PMF we have assumed that planets in the HZ of late K and dM stars (M < 0.7M ⊙ ) are tidally locked at times t < 0.7 Gyr (n=2 in eq. (33), see figure 5). Planets around G and early K stars (M 0.7M ⊙ ) will be assumed to have their primordial periods of rotation that we chose in the range 1 − 100 days as predicted by models of planetary formation (Miguel & Brunini 2010). Figures 7 and 8 show the evolution of magnetosphere properties for tidally locked and unlocked potentially habitable planets respectively. In all cases we have assumed that the planets are in the middle of the HZ of their host stars.
Even at early times tidally locked planets of arbitrary mass have a non-negligible magnetosphere radius R S > 1.5 R p . Previous estimates of the standoff distances for tidally locked planets are much lower than the values reported here. As an example Khodachenko et al. 2007 place the standoff distances well below 2 R p , even under mild stellar wind conditions (see figure 4 in their work) and independent of planetary mass and age. In contrast our model predicts standoff distances for tidally locked planets in the range of 2-6 R p depending on planetary mass and stellar age. The differences between both predictions arise mainly from the underestimation of the dipole moment for slowly rotating planets found in these works. Thermal evolution and the dependency on planetary mass of the PMF properties are responsible for the rest of the discrepancies in previous estimates of the magnetosphere properties.
Though tidally locked planets seem to have larger magnetospheres than previously expected, they still have large polar caps, a feature that was previously overlooked. As a consequence of this fact well protected atmospheres, i.e. atmospheres that lie well inside of the magnetosphere cavity (hereafter magnetised planets), could have more than 15% of their surface area exposed to open field lines where thermal and non-thermal processes could efficiently remove atmospheric gases. Moreover, our model predicts that these planets would have multipolar PMF which contributes to an increase of the atmospheric area open to the interplanetary and magnetotail regions (Siscoe & Crooker 1976). Then the exposition of magnetised planets to harmful external effects would be a complex function of the standoff distance and the polar cap area.
Overall magnetic protection improves with time. As the star evolves the dynamic pressure of the stellar wind decreases more rapidly than the dipole moment (see figures 4 and 6). As a consequence the standoff distance grows in time and the polar cap is shrunk. However with the reduction in time of the stellar wind pressure the magnetopause field is also reduced a fact that can affect the incoming flux of CR at late times.
The sinuous shape of the contour lines in the middle and lower rows is a byproduct of the inner core solidification in planets with M p < 2M ⊕ . Critical boundaries between regions with very different behaviors in the magnetosphere properties are observed at M p ∼ 1.0M ⊕ and M p ∼ 1.8M ⊕ in the middle and rightmost panels of the standoff radius and polar cap area contours. Planets to the right of these boundaries still have a completely liquid core and therefore produce weaker PMFs (lower standoff radius and larger polar cap areas). On the other hand, the inner core in planets to the left of the these boundaries have already started to grow and therefore their PMFs are stronger.
Unlocked planets (figure 8) are better protected than slowly rotating tidally locked planets by developing extended magnetospheres R S 4 R p and lower polar cap areas A pc 10%. It is interesting to notice that in both cases and at times t ∼ 1 Gyr a smaller planetary mass implies a lower level of magnetic protection (lower standoff distances and larger polar caps). This result seems to contradict the idea that low-mass planets (M p 2) are better suited to develop intense and protective PMFs (Gaidos et al. 2010;Tachinami et al. 2011;Zuluaga & Cuartas 2012). To explain this contradiction one should take into account that magnetic protection as defined in this work depends on dipole moment instead of surface magnetic field strength. Since dipole moment scales-up as M ∼ B dip R 3 p more massive planets will have a better chance to have large and protective dipole moments.
It is interesting to compare the predicted values of the maximum dipole moment calculated here with the values roughly estimated in previous attempts (Grießmeier et al. 2005 In their, work López-Morales et al. 2011 use the same power-based scaling laws we applied here but assume a rather simple interior model and a static thermal model where the convective power is set such that maximizes the efficiency with which the convective energy is converted into the magnetic field. A more detailed account of the evolution of magnetosphere properties for the already discovered habitable planets is presented in figure 9. In all cases we have assumed that all planets have compositions similar to Earth (RTEM). Although almost all planets are tidally locked, we have also computed the magnetic properties for a primordial period of rotation P = 1 day.
The case of the "hydrated" Venus is particularly interesting in order to analyse the rest of planets. The dynamo of the actual Venus probably shut down at t = 3 Gyr as a consequence of the drying of the mantle (Christensen et al. 2009)   induced by a runaway greenhouse and insufficient early magnetic protection played a central role in the extinction of the early Venusian PMF. The evolution of the PMF in the potentially habitable planets GL 581d, GJ 667Cc and HD 40307g could have a similar fate. Their masses are much larger and therefore their atmospheres are protected by stronger gravitational fields. For planet HD 40307g our reference thermal evolution model predicts a late shut down of the dynamo t dyn ∼ 4 Gyr. According to our reference model the planet is presently devoid of a dynamo generated magnetic field. However, being around a K star (M ⋆ ∼ 0.7) the stellar wind and XUV radiation have probably decreased enough to not represent at present times a real threat for its atmosphere.

. A massive loss of water
GL 581d and GJ 667Cc are located in the HZ of dM stars where the stellar wind pressure and XUV radiation, even at times as late as 4 Gyr, are intense enough to erode their atmospheres or to make them lose their water content. The RTEM predicts that for an Earth-composition GL 581d at present times had already lost its dynamo and has been exposed for almost 2.5 Gyr to the harmful effects of the stellar wind and CR. The planet however is the most massive of the three planets and probably has a thick atmosphere able to withstand the continuous agression of its host star.
Given the estimated age of the GJ 667C system (t ≈ 2 Gyr), the RTEM predicts that the planet still has a dynamo (red circle in figure 9). Magnetosphere properties are very close to that of our "hydrated" Venus, 4 Gyr ago. However its mass is lower than that of GL 581d and it is located at the inner edge of the HZ where the exposition to the XUV radiation from its host star (a young M1 star) could have been enough to induce massive loss of atmospheric gases including water. We will come back on these issue in section 6.1 when we will show how to estimate the atmospheric mass-loss rate for this particular planet.

Toward an estimation of the atmospheric mass-loss
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 atmopsheres of magnetised and unmagnetised potentially habitable planets. This is a fundamental goal to be pursued in the near future if we want to assess the actual habitability of present and future discovered terrestrial planets in the HZ of their host stars. 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 in the 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 2009and 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 2 rich atmospheres (Watson et al. 1981;Kulikov et al. 2006;Tian et al. 2008) and dry Venus-like CO 2 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 nonthermal 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 or larger than the magnetic standoff distance R S we will say that the planet is unmagnetised. 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) atmopsheric 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 unmagnetised 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 = 35 − 10 PEV (1 PEV = 0.64 erg cm −2 s −1 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 PEV while GJ 667Cc received the maximum amount of XUV radiation among them, F XUV = 450 − 800 PEV.
Using the recent results by Tian 2009 that computed the exosphere properties of massive super-Earths, i.e. M p 6M ⊕ , subject to different XUV fluxes, we can estimate the exosphere radius and mass-loss rate for our three habitable super-Earths. Actually, since the Tian (2009) results are only available for planets with a minimum mass of M p = 6M ⊕ , a qualitative extrapolation of the results for 10, 7 and 6 M ⊕ (see figure 4 and 6 in its paper), shows that exobase radius and mass-loss rates are larger for less massive planets. This is particularly useful at trying to apply the Tian's results to GJ 667Cc M p ≈ 4.5M ⊕ 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 ≈ 7M ⊕ ) and assuming an initial CO 2 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 predict the planet was enshrouded by a protective magnetosphere (see figure 9). After dynamo shut down 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 been eroding the atmosphere for 3-4 Gyrs. During this unmagnetised phase the atmospheric mass-loss rate can be simply estimated asṀ ≈ αmnv eff (Zendejas et al. 2010) where α is the socalled 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 eq. (2)). Using a entrainment efficiency α ∼ 0.3 (which is appropriate for example to describe the mass-loss of the Venus atmosphere) we obtain that the total mass-loss during the unmagnetised 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 preserve a dense enough atmosphere able to sustain surface liquid water and hence to be actually habitable.
The case of GL 581d (M p ≈ 6M ⊕ ) is quite different. Assuming that our estimations of the XUV flux are right, 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 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 atmopsheric composition and of course in the magnetic model developed here should lead to a different conclusion and further theoretical and probably observational investigation is required.
The most interesting case is that of GJ 667Cc (M p ≈ 4.5M ⊕ ). The minimum exosphere radius predicted for this planet at t ∼ 1 Gyr lies between 3.0 − 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 2 rich atmosphere will be larger than 2 − 4 × 10 10 atoms cm −2 s −1 (see figure 6 in Tian 2009). We should recall that this is actually the value for a 6M ⊕ super-Earth. For the actual mass of the planet, 4.5M ⊕ , the mass-loss rate could be even larger. Moreover if as predicted here the exosphere is exposed directly to the stellar wind, non-thermal processes can contribute to a larger increase in the mass-loss from the planetary atmosphere. At the minimum massloss rate the exposed GJ 667Cc atmosphere could have lost more than ∼ 10 46 atoms of Carbon in just ∼ 100 Myr and in the first gigayear the amount of carbon thermally lost to space could rise to ∼ 10 47 atoms. If we scale-up linearly with planetary mass the total inventory of CO 2 in the atmosphere, crust and mantle of Venus, which is 2 − 3 × 10 46 molecules (see Tian 2009 and references therein), a 4.5M ⊕ planet will have a total budget of ∼ 10 47 CO 2 molecules. In summary at the minimum mass-loss rate and assuming a relatively rapid degassing of the planet, Gj 667C c could have lost its total inventory of Carbon to interplanetary space in the first couple of gigayears. Even assuming that large amounts of CO 2 are still trapped in the mantle and crust of the planet, its atmosphere should be being rapidly obliterated by the stellar wind. We speculate that GJ 667Cc is a sort of "Venus-like" planet. Regardless the fact the planet is well inside the radiative habitable zone it has lost its capacity to support life via a massive stellar-wind induced loss of volatiles.

DISCUSSION
Applying simplified thermal evolution model and dynamo scaling laws to planets whose bulk properties are barely known or even hypothetical it is challenging and probably raises more questions than it attempts to respond. Further observations of the potentially habitable planets should be required to precise their physical properties and to reliably model its interior and thermal evolution. Moreover continued observational efforts to look for direct evidence or proxies of planetary magnetospheres and any other signatures of magnetic protection, though challenging, should also be attempted. Here we want to discuss the assumptions on which our global model relies, its uncertainties as measured by the sensitivity of the model to changes in key physical parameters and the missing pieces of information and present observational limitations to confirm or improve this and other models of planetary magnetic protection.

Model assumptions
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 we discuss in the following paragraphs.
We have assumed that terrestrial planets 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 bestcase 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 eq. (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 eq. (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 eq. (3), is equivalent to substracting it from the stellar wind dynamic pressure P sw in the right-hand side. An effective stellar wind pressure P ′ sw = P sw (1 − P pl /P sw ) would replace the stellar wind term in the standoff distance definition (eq. (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 be increased by only a 10%. On the other hand in order to have a standoff distance one order of magnitude larger than that estimated using eq. (5) the plasma pressure inside the magnetosphere should amount for 99.999% of the total pressure. This is precisely what an unmagnetised planet would look like. In summary including more realistic condition into the definition of the magnetosphere boundary will not modify too much our results.
A final but not 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 (CME). To model the effect that a steady flux/influx of CME plasma could have in closein 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 6 . Taking into account that R S ∼ P −1/6 SW we found that under the harsher conditions the magnetosphere radius and polar cap areas will be modified only by 10-30% in 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.

Sensitivity analysis
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 6 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 already discovered habitable planets (see table 2) including a hydrated Venus and present Earth. In all cases we assume for simplicity a primordial period of rotation of P = 1 day.
Since the dipole moment is an evolving quantity we have plotted in figure 10 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 has 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 6 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 have 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 the Earth's value CMF = 0.325, i.e. assumes that all planets are dominated by a sillicate-rich mantle. As a comparison Mars has a CMF = 0.23 and the value for Mercury is CMF = 0.65 (it should be recalled that Mercury could have lost a significant fraction of its mantle sillicates 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 eq. (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 sensitivty 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 provides at least the magnitude and sign of the effect that this parameter has in the dynamo properties predicted by our thermal evolution model.

3.
High pressure viscosity rate coefficient, b (see eq. (26)). Rheological properties of sillicates 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 eqs. (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 eq. (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. In this value all the figures are significative reflecting 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 out from the core. In the RTEM we used a value k c =40 Wm −1 K −1 that fits the thermal evolution history and present magnetic field of the Earth (see table 1). Although recent first-principles analysis suggest that values as large as 150-250 Wm −1 K −1 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 that found by Pozzo et al. 2012 should be attempted. For our sensitivity analysis we varied k c between 35 and 70 Wm −1 K −1 , 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 specially at core conditions. It affects strongly 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 kind of core alloys a relatively large range of values of this parameter has been used in literature (see VAL06 and references there in). 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 , that is also varied to study the sensitivity of thermal evolution models , were also studied. We did not found significant sensitivity of our model to variation of those parameters.
We depict in figure 10 the relative variation in the aforementioned magnetosphere and dynamo properties when each of the previously described parameters were varied independently.
We have found that planetary composition (CMF), 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 core surface is of the same order for Fe-poor and Fe-rich planets. However a small core means also a lower magnetic dipole moment, i.e. M ∼ R 3 c . Planets with lower content of iron also have small and hot cores and therefore the solid inner core formation and the shutting down of 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 earth-like planets. A low viscosity lower mantle will favour the extraction of heat from the metallic core increasing the available convective energy for dynamo action. On the other hand a viscuous 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 are negligible, at least on what 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 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 with a larger Grüneisen parameter hand, earthlike planets could get a solid inner-core very early in their thermal histories even almost since the beginning. This is the 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 produce 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 < 2M ⊕ 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 a colder core develops a solid inner core almost from the beginning and the release of latent and gravitational energy feeds a stronger dynamo.
The case of more massive planets, M p > 2M ⊕ , 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 affect less the results of the thermal evolution model. Besides the case of massive planets where differences in the order of 10-30% in the magnetic properties are observed when k c varies between its extremes, the magnetic properties calculated with our reference thermal evolution model seem very robust against variations in 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 in the order of magnitude. Moreover since the standoff distance and polar cap areas, which are the actual proxies to magnetic protection, goes as M 1/3 , a one order of magnitude estimation of M will give us an estimation of the level of magnetic protection off by a factor around 2.
To clarify this point let's 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 notice that probably 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.

Observational support
Validating or improving thermal evolution models and dynamo scaling laws for the case of super-Earths nowadays represents a huge observational challenge. The available sensitivity of our best instruments in the Earth and in space are rather insufficent. New and/or improved instruments and observational techniques will be required to detect, catalogue and compare the ther-mal and magnetic properties of low mass planets in the medium to far future. However the importance that the detection and characterization of the magnetic properties of future discovered potentially habitable planets in order to assess their true habitability clearly justify the effort.
The first goal seems to be the direct or indirect detection of super-Earth magnetospheres. Four methods, some of them 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 (radioemission) 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. 2010Vidotto et al. , 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 terrestrial planets 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, LOFAR or the Long Wavelength Array LWA, the expected power and spectra are several orders of magnitude below the threshold of detection. Powers as large as 10 3 -10 5 times the Jupiter radioemission 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 super-Earths magnetospheres are several orders of magnitude lower than these thresholds and probably are far from being 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. Vidoto (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 should be required to have a detectable signal of a bow shock (A.A. Vidotto, Private Communication). Our habitable super-Earths have magnetopause fields in the order of a few micro Gauss (see figure 7). The case of an ion tail coming from a weakly magnetised planet has received less attention and probably could offer better chances for a future indi-rect 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 study directly and indirectly planetary magnetospheres. Planetary auroras with intensities as high as 10 2 to 10 3 times larger than that of the Earth are expected in closein 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 −11 photons m −2 s −1 srad −1 ) 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 13 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 habitable zone, the total FUV power produced by auroras in the polar cap of earth-like planets could be only 10 5 W which is 8 orders of magnitude less than the present detection threshold. Not to mention that the FUV radiation should be detected against an intense UV background coming from a probably young and active low mass star. If we can find ways to overpass these difficulties, the observation of the FUV and X-ray emmision and its variability from auroras in potentially habitable super-Earths could be used as powerful probes of the magnetic environment around the planet.

SUMMARY AND CONCLUSIONS
We studied here the influence that the thermal evolution of potentially habitable terrestrial planets has in the protection that an evolving planetary magnetopshere could provide against the atmospheric erosion caused by the stellar wind.
We developed a simple parametrized thermal evolution model able to reproduce the global thermal history and magnetic properties of the past and present Earth. We applied this model to predict the thermal histories of planets with masses in the range of 0.5 to 6.0 M ⊕ and with chemical compositions similar to Earth. Using these results and applying up-to-date dynamo scaling laws we predicted the magnetic properties of terrestrial planets in the habitable zone as a function of time, planetary mass and rotation rate. A simple model of the evolution and interaction of the stellar wind with the planetary magnetic field, that has been adapted from previous works, allowed us to compute the global properties of the magnetosphere in order to assess the level of magnetic protection that potentially habitable Earth-like planets could actually have. We applied our model to the case of already known potentially habitable terrestrial planets (GL 581d, HD 40307g and GJ 667Cc), to the Earth itself and to the case of a hydrated Venus. In the case of the Earth our model reproduces fairly well the early and present thermal and magnetic properties of our planet. In the case of the hydrated Venus, the model predicts low values the standoff distance and large polar cap areas in the frist critical gigayear of planetary and thermal evolution, which are compatible with the idea that the planet lacked a strong enough magnetic protection able to avoid a massive loss of water and volatiles that finally lead to the shut down of its dynamo ∼ 3 Gyr ago.
Compelling results were found in the case of the three already discovered extra-solar-system potentially habitable planets. Assuming an earth-like composition and thermal evolution parameters similar to those used in the case of the Earth (reference thermal evolution model, RTEM), our model predicted that the dynamo of GL 581d and HD 40307g have been already shut down. A younger GJ 667Cc seems to still have an active dynamo.
A non-trivial dependence of the magnetic properties on planetary age, planetary mass and period of rotation has been found in general for terretrial planets inside the HZ of their host stars. Thermal evolution is responsible for the non-trivial relationship among all these properties. Contrary to what was found in previous work, tidally locked planets could develop relatively intense magnetic fields and extended magnetospheres. However they also have extended polar caps and probably multipolar magnetic fields where field lines open to the interplanetary space and magnetotail regions probably increasing the non-thermal mass-losses.
Using recent results for the relationship between the exposition to XUV radiation, the exobase radius and mass-loss rate from massive super-Earths, we estimated the level of exposure and mass-losses for the three already discovered potentially habitable super-Earths. With the available information not too much could be said about the magnetic protection of HD 40307g. Further theoretical investigations are required to evaluate this case. Our model predicts a large enough magnetosphere able to protect GL 581d against the erosive action of the stellar wind during the first critical phases of planetary evolution. However since our model is still optimistic further theoretical and probably observational analyses should be performed to establish on a more solid basis the magnetic protection of this planet. Our upper-limit to the standoff-distance and the most optimistic estimation of exobase radius and mass-loss rate from the atmosphere of GJ 667Cc, point-out the fact that this planet has already lost a large fraction of its inventory of volatiles. All the evidence compiled in this work make GJ 667Cc a sort of "Venus analogue". Although further theoretical analyses are required our best guess is that despite the fact that it is inside the radiative HZ of its host star the planet is presently uninhabitble.
Last but not least we tested the robustness of our conclussions by changing several of the most sensitive input parameters of our thermal evolution model. We found that even under the present uncertainties the predicted properties of planetary magnetopsheres are rather robust. We calculated that introducing large variations in the composition of the planets and the rheological and thermal properties of their interiors with respect to the reference thermal evolution model, the critical magnetic properties, such as the standoff radius and the area of the polar cap, change only by a factor of two. Results are also robust against uncertainties in the stellar wind properties that could be very important in the case of close-in habitable planets around active and young dM stars.
The problem of evaluating the magnetic protection of potentially habitable planets is far from being settled. Other sources of intrinsic magnetic fields, thermal evolution and interior structure of planets with "exotic" compositions, improved theoretical models and new experi-mental evidence of the behavior of iron at high pressures and temperatures, improved and validated models of the evolution and spatial structure of stellar winds and of course more and better observational data coming from the already discovered habitable super-Earths and future discovered potentially habitable exoplanets, will allow us to assess the actual magnetic protection of potentially habitable environments.
We appreciate the useful discussion and comments of Mercedes Lopez-Morales and other colleagues participating in the Exploring Strange New Worlds 2011 Conference (Arizona, U.S.) and in the VI Taller de Ciencias Planetarias 2012 (Montevideo, Uruguay). Special thanks to Ignacio Ferrin who with his clever questions and suggestions originally motivate us to pursue some of the goals of this work. We also thank to Lisa Kaltenegger and Jeffrey Linsky for their insightful comments on preliminary versions of this manuscript. We want to give special thanks to all our fellow colleagues abroad that have provided us with some key literature unobtainable from our country. We are grateful to Aaron West and Luke Webb for his careful revision of the English in the manuscript. The remnant errors are all ours. Anonymous referee contributed significantly not only to the improvement of the manuscript but to the quality of the research conclussions and should be also aknowledged. PC is supported by the Vicerrectoria de Docencia of the University of Antioquia. This work has been completed with the financial support of the CODI-UdeA under the grant IN591CE and by the University of Medellin under the grant number 530.   (18)). Middle panel: radius of the inner core R ic . Lower panel: time of inner core formation (blue squares) and dynamo lifetime (red circles). In the RTEM the metallic core is liquid at t = 0 for all planetary masses. Planets with a mass Mp < M crit = 2.0M ⊕ develop a solid inner core before the shut down of the dynamo while the core of more massive planets remains liquid at least until the dynamo shutdown.  (29) and (31) for TPs 0.5-4.0 M ⊕ . We plot the local Rossby number (lower panel), the maximum surface dipole field (middle panel) and the maximum dipole moment (upper panel). We included the present values of the geodynamo (⊕ symbol) and three measurements of paleomagnetic intensities (error bars) at 3.2 and 3.4 Gyr ago (Tarduno et al. 2010). We compare the magnetic properties for two periods of rotation, 1 day (solid curves) and 2 days (dashed curves). The effect of a larger period of rotation is more significant at early times in the case of massive planets (Mp 2M ⊕ ) and at late times for lower mass planets.  (Zuluaga & Cuartas 2012): rapid rotating planets (P 1 day), dipole moments are large and almost independent of rotation rate; slowly rotating planets (1 day P 5 day), dipole moments are intermediate in value and highly dependent on rotation rate; and very slowly rotating planets (P 5 − 10 days), small but non-negligible rotation-independent dipole moments. For (Mp < 2M ⊕ ) the shape of the dipole-moment contours is determined by t ic .  Kopparapu et al. (2013). Stellar properties are computed at τ = 3 Gyr using the models by Baraffe et al. 1998. Planets at distances below the dashed line would be tidally locked before 0.7 Gyr (Peale 1977). The location of Earth, Venus and the potentially habitable extra-solar-system planets GL 581d, GJ 667 C c and HD 40307g are also included. Fig. 6.-Evolution of the stellar wind dynamic pressure at the center of the HZ for a selected set of stellar masses. The reference average solar wind pressure is P SW⊙ = 1.86 nPa. Dashed curves indicate the value of the stellar wind pressure at the inner and outer edges of the HZ around stars with 0.2 M ⊙ and 1.0 M ⊙ respectively. The HZ limits where the pressure were calculated are assumed static and equal to those at τ = 3 Gyr. Stellar wind pressure at t < 0.7 Gyr computed with the semiempirical model used in this work is too uncertain and were not plotted. Fig. 7.-Evolution of the magnetopause field (upper row), standoff distance (middle row) and polar cap area (lower row) of tidally locked (slow rotating) planets around late dK and dM stars. The rotation of each planet is assumed equal to the orbital period at the middle of the HZ (see values in the rightmost vertical axes). The value of the magnetosphere properties rerturned by the contour lines in these plots could be an under or an overestimation of these properties according to the position of a planet inside the HZ. In the case of GL 581d (GJ 667Cc), which is located in the outer (inner) edge of the HZ, the magnetopause field and polar cap area are overestimated (underestimated) while the standoff distance is underestimated (overestimated). . For all planets we have assumed a constant period of rotation P = 1 day. Fig. 9.-Evolution of magnetosphere properties for the already discovered habitable SEs, the Earth and an "hydrated" version of Venus (low viscosity mantle and mobile lid). Shaded regions are bounded by the properties calculated at a minimum period of rotation of P ≈ 1 day (upper and lower bounds in standoff radius and polar cap area curves respectively) and the maximum period of rotation P ≈ Po corresponding to a perfect match between the rotation and orbital periods (tidal locking). Magnetopause fields do not depend on the rotation period of the planet. Filled circles are the predicted present day magnetosphere properties computed according to the properties summarized in table 2. Fig. 10.-Sensitivity analysis of our reference thermal evolution model (RTEM). Squares and diamonds indicate the relative value of three critical magnetic properties, M dip (the average of the dipole moment between 0.7 and 2 Gyr), t ic (time of inner core formation) and t dyn (dynamo lifetime), as calculated by the thermal evolution model. For the analysis 5 different key thermal evolution parameters were independently changed with respect to the reference value in the RTEM: the core mass fraction (CMF), the thermal conductivity of the core (kc), the Grüneisen parameter at core conditions (γ 0c ), the high pressure viscosity rate coefficient (b) and the adiabatic factor (ǫ adb ). The results obtained when the minimum value of the parameters were used are indicated with squares. Conversely, the results obtained with the maximum value of each parameter are indicated with diamonds.