Conservative motion

The Hamiltonian of the motion can be written using the classical canonical Andoyer's action and angle variables (Andoyer Reference Andoyer1923; Kinoshita Reference Kinoshita1977), but these variables become degenerate due to a simplifying assumption that is done here: the planet is considered to be a rigid body with moments of inertia A=B<C, and we merge the figure axis with the direction of the angular momentum (gyroscopic approximation)Footnote ¹. We are then left with only two action variables (L, X) and their conjugate angles (θ, φ) (Néron de Surgy & Laskar Reference Néron de Surgy and Laskar1997)Footnote ². The quantity L=Cω, with rotation rate ω, is the the angular momentum along the C-axis, and X=L cos ε, where ε is the obliquity, is its projection on the normal to the quasi-inertial ecliptic plane. The variable θ is the hour angle between the equinox and a fixed point of the equator, and φ is the precession angle (Fig. 1). Averaging the Hamiltonian over the rotation angle and the mean anomaly, we getFootnote ³ (Kinoshita Reference Kinoshita1977; Néron de Surgy & Laskar Reference Néron de Surgy and Laskar1997; Correia & Laskar Reference Correia and Laskar2010b)

(1)$$\bar H = \displaystyle{{L^2} \over {2C}} - {\rm \alpha} \displaystyle{{X^2} \over {2L}}\;, $$

where

(2)$${\rm \alpha} = \displaystyle{{3GM_*} \over {2a^3 (1 - e^2 )^{3/2}}} \displaystyle{{E_{\rm d}} \over {\rm \omega}} \simeq \displaystyle{3 \over 2}\displaystyle{{n^2} \over {\rm \omega}} \left( {1 + \displaystyle{3 \over 2}e^2} \right)E_{\rm d} $$

is the ‘precession constant’, while quantities G, M _*, a, n and e are the gravitational constant, stellar mass, semi-major axis, mean motion and eccentricity, respectively. Quantity E _d is the dynamical ellipticityFootnote ⁴,

(3)$$E_{\rm d} = \displaystyle{{C - A} \over C} = \displaystyle{{k_{\rm f} R^5} \over {3GC}}{\rm \omega} ^2 + {\rm \delta} E_{\rm d}, $$

where R is the planet radius, and k _f is the fluid Love number. The first part of this expression corresponds to the flattening in hydrostatic equilibrium (Lambeck Reference Lambeck1980), and the second corresponds to the departure from this equilibrium.

Fig. 1. Andoyer's canonical variables. L is the projection of the total rotational angular momentum vector L on the principal axis of inertia k, and X the projection of the angular momentum vector on the normal to the orbit K. The angle between the line of nodes γ and a fixed point of the equator A is the hour angle, θ, whereas the angle between a reference point in the orbit and the line of nodes γ is the precession angle, φ.

Since Andoyer's variables are canonical, the spin equations of motion are easily obtained from the mean Hamiltonian (equation (1)) as

(4)$$\displaystyle{{{\rm d}L} \over {{\rm d}t}} = - \displaystyle{{{\rm \partial} \bar H} \over {\partial {\rm \theta}}} \;, \quad \displaystyle{{{\rm d}X} \over {{\rm d}t}} = - \displaystyle{{\partial \bar H} \over {\partial {\rm \varphi}}} \;, \quad \displaystyle{{{\rm d}{\rm \varphi}} \over {{\rm d}t}} = \displaystyle{{\partial \bar H} \over {\partial X}},$$

which gives

(5)$$\displaystyle{{{\rm d}L} \over {{\rm d}t}} = \displaystyle{{{\rm d}X} \over {{\rm d}t}} = 0,\quad {\rm and}\quad \displaystyle{{{\rm d}{\rm \varphi}} \over {{\rm d}t}} = - {\rm \alpha} \cos {\rm \varepsilon}, $$

i.e., the rotation rate and the obliquity are constant, and the planet precesses at a constant rate.

Tidal effects

Tidal effects arise from planetary differential and inelastic deformations caused by a perturbing body. There are two types of tidal effects: the gravitational tides and the thermal atmospheric tides. The estimations for both effects are based on a general formulation of the tidal potential, initiated by Darwin (Reference Darwin1880).

Gravitational tides

Gravitational tides are raised on the planet by a perturbing body because of the gravitational gradient across the planet. The force experienced by the side facing the perturbing body is stronger than that experienced by the far side. These tides are mainly important upon the solid (or liquid) part of the planet, and are independent of the existence of an atmosphere.

Since the planets are not perfectly rigid, there is a distortion that gives rise to a tidal bulge. This redistribution of mass modifies the gravitational potential generated by the planet in any point of the space. The additional amount of potential, the tidal potential U _g, is responsible for the modifications in the planet's spin (and orbit), and it is given byFootnote ⁵ (e.g., Lambeck Reference Lambeck1980):

(6)$$U_{\rm g} = - k_2 \displaystyle{{GM_*^2} \over R}\left( {\displaystyle{R \over {r_*}}} \right)^3 \left( {\displaystyle{R \over r}} \right)^3 P_2 (\cos \,S),$$

where r and r _* are the distances from the planet's centre of mass to a generic point and to the star, respectively, S is the angle between these two directions, P ₂ are the second-order Legendre polynomials, and k ₂ is the second potential Love number.

Expressing the tidal potential given by expression (6) in terms of Andoyer angles (θ, φ), we can obtain the contribution to the spin evolution from expressions (4) using U _g at the place of ${\overline {\cal H}}$. As we are interested here in the study of the secular evolution of the spin, we also average U _g over the periods of mean anomaly and longitude of the periapse of the orbit. This work is done with the help of the algebraic manipulator TRIP (Laskar Reference Laskar1989, Reference Laskar1994), which expands the potential in Fourier series, as in Kaula (Reference Kaula1964) and Correia & Laskar (Reference Correia and Laskar2010b). For a planet orbiting its host star, where the star is both the perturbing and interacting body (r=r _*), we then find for the averaged equations of motion:

(7)$$\displaystyle{{{\rm d}L} \over {{\rm d}t}} = K_{\rm g} \mathop \sum \limits_{\rm \sigma} b_{\rm g} ({\rm \sigma} ){\rm \Lambda} _{\rm \sigma} ^{\rm g} (x,e),$$

(8)$$\displaystyle{{{\rm d}X} \over {{\rm d}t}} = K_{\rm g} \mathop \sum \limits_{\rm \sigma} b_{\rm g} ({\rm \sigma} ){\rm \Upsilon} _{\rm \sigma} ^{\rm g} (x,e),$$

where

(9)$$K_{\rm g} = - \displaystyle{{GM_*^2 R^5} \over {a^6}}, $$

and x=X/L=cos ε. The coefficients, Λ_σ^g and ${\rm \Upsilon} _{\rm \sigma} ^{\rm g} $, are polynomials in the eccentricity (Kaula Reference Kaula1964). For a planet with moderate eccentricity, we may neglect terms in e ⁴ and greater, and obtain the following expansion for the previous equations:

(10)$$\eqalign{\displaystyle{1 \over {K_{\rm g}}} \displaystyle{{{\rm d}L} \over {{\rm d}t}} = & b_{\rm g} (2{\rm \omega} + 3n)\displaystyle{{147} \over {128}}e^2 (1 - x)^4 \cr & + b_{\rm g} (2{\rm \omega} + 2n)\displaystyle{3 \over {32}}(1 - 5e^2 )(1 - x)^4 \cr & + b_{\rm g} (2{\rm \omega} + n)\displaystyle{3 \over {128}}e^2 (37 + 70x + 37x^2 )(1 - x)^2 \cr & + b_{\rm g} (2{\rm \omega} )\displaystyle{3 \over 8}(1 + 3e^2 )(1 - x^2 )^2 \cr & + b_{\rm g} (2{\rm \omega} - n)\displaystyle{3 \over {128}}e^2 (37 - 70x + 37x^2 )(1 + x)^2 \cr & + b_{\rm g} (2{\rm \omega} - 2n)\displaystyle{3 \over {32}}(1 - 5e^2 )(1 + x)^4 \cr & + b_{\rm g} (2{\rm \omega} - 3n)\displaystyle{{147} \over {128}}e^2 (1 + x)^4 \cr & + b_{\rm g} ({\rm \omega} + 3n)\displaystyle{{147} \over {64}}e^2 (1 - x^2 )(1 - x)^2 \cr & + b_{\rm g} ({\rm \omega} + 2n)\displaystyle{3 \over {16}}(1 - 5e^2 )(1 - x^2 )(1 - x)^2 \cr & + b_{\rm g} ({\rm \omega} + n)\displaystyle{3 \over {64}}e^2 (1 - 2x + 37x^2 )(1 - x^2 ) \cr & + b_{\rm g} ({\rm \omega} )\displaystyle{3 \over 4}(1 + 3e^2 )(1 - x^2 )x^2 \cr & + b_{\rm g} ({\rm \omega} - n)\displaystyle{3 \over {64}}e^2 (1 + 2x + 37x^2 )(1 - x^2 ) \cr & + b_{\rm g} ({\rm \omega} - 2n)\displaystyle{3 \over {16}}(1 - 5e^2 )(1 - x^2 )(1 + x)^2 \cr & + b_{\rm g} ({\rm \omega} - 3n)\displaystyle{{147} \over {64}}e^2 (1 - x^2 )(1 + x)^2,} $$

and

(11)$$\eqalign{\displaystyle{1 \over {K_{\rm g}}} \displaystyle{{{\rm d}X} \over {{\rm d}t}} = & - b_{\rm g} (2{\rm \omega} + 3n)\displaystyle{{147} \over {128}}e^2 (1 - x)^4 \cr & - b_{\rm g} (2{\rm \omega} + 2n)\displaystyle{3 \over {32}}(1 - 5e^2 )(1 - x)^4 \cr & - b_{\rm g} (2{\rm \omega} + n)\displaystyle{3 \over {128}}e^2 (1 - x)^4 \cr & + b_{\rm g} (2{\rm \omega} - n)\displaystyle{3 \over {128}}e^2 (1 + x)^4 \cr & + b_{\rm g} (2{\rm \omega} - 2n)\displaystyle{3 \over {32}}(1 - 5e^2 )(1 + x)^4 \cr & + b_{\rm g} (2{\rm \omega} - 3n)\displaystyle{{147} \over {128}}e^2 (1 + x)^4 \cr & - b_{\rm g} ({\rm \omega} + 3n)\displaystyle{{147} \over {32}}e^2 (1 - x^2 )(1 - x)^2 \cr & - b_{\rm g} ({\rm \omega} + 2n)\displaystyle{3 \over 8}(1 - 5e^2 )(1 - x^2 )(1 - x)^2 \cr & - b_{\rm g} ({\rm \omega} + n)\displaystyle{3 \over {32}}e^2 (1 - x^2 )(1 - x)^2 \cr & + b_{\rm g} ({\rm \omega} - n)\displaystyle{3 \over {32}}e^2 (1 - x^2 )(1 + x)^2 \cr & + b_{\rm g} ({\rm \omega} - 2n)\displaystyle{3 \over 8}(1 - 5e^2 )(1 - x^2 )(1 + x)^2 \cr & + b_{\rm g} ({\rm \omega} - 3n)\displaystyle{{147} \over {32}}e^2 (1 - x^2 )(1 + x)^2 \cr & - b_{\rm g} (3n)\displaystyle{{441} \over {64}}e^2 (1 - x^2 )^2 \cr & - b_{\rm g} (2n)\displaystyle{9 \over {16}}(1 - 5e^2 )(1 - x^2 )^2 \cr & - b_{\rm g} (n)\displaystyle{9 \over {64}}e^2 (1 - x^2 )^2.} $$

The coefficients b _g(σ) are related to the dissipation of the mechanical energy of tides in the planet's interior, responsible for a time delay Δt _g(σ) between the position of ‘maximal tide’ and the sub-stellar point. They are related to the geometric lag δ_g(σ) as:

(12)$$b_{\rm g} ({\rm \sigma} ) = k_2 \sin {\rm \;} 2{\rm \delta} _{\rm g} ({\rm \sigma} ) = k_2 \sin ({\rm \sigma} {\rm \Delta} t_{\rm g} ({\rm \sigma} )).$$

Dissipation equations (10) and (11) must be invariant under the transformation (ω, x) by (−ω, −x) which implies that b(σ)=−b(−σ). That is, b(σ) is an odd function of σ (Correia et al. Reference Correia, Laskar and Néron de Surgy2003). Although mathematically equivalent, the couples (ω, x) and (−ω,−x) correspond to two different physical situations (see Correia & Laskar Reference Correia and Laskar2001).

Thermal atmospheric tides

The differential absorption of the Solar heat by the planet's atmosphere gives rise to local variations of temperature and consequently to pressure gradients. The mass of the atmosphere is then permanently redistributed, adjusting for an equilibrium position. More precisely, the particles of the atmosphere move from the high-temperature zone (at the sub-solar point) to the low-temperature areas (e.g., Arras & Socrates Reference Arras and Socrates2010). Observations on Earth show that the pressure redistribution is essentially a superposition of two pressure waves: a diurnal tide of small amplitude and a strong semi-diurnal tide (Bartels Reference Bartels1932; Haurwitz Reference Haurwitz1964; Chapman & Lindzen Reference Chapman and Lindzen1970).

As for gravitational tides, the redistribution of mass in the atmosphere gives rise to an atmospheric bulge that modifies the gravitational potential generated by the atmosphere in any point of the space. The tidal potential U _a responsible for the spin changes is given byFootnote ⁶ (e.g., Correia & Laskar Reference Correia and Laskar2003a):

(13)$$U_{\rm a} = - \displaystyle{3 \over 5}\displaystyle{{\tilde p_2} \over {{\rm \bar \rho}}} \left( {\displaystyle{R \over r}} \right)^3 P_2 (\cos \,S),$$

where $\tilde p_2 $ is the second-order surface pressure variations, and ${\rm \bar \rho} $ is the mean density of the planet.

To find the contributions to the spin we use expressions (4) together with the averaging method over fast varying angles (mean anomaly and longitude of the periapse), which gives:

(14)$$\displaystyle{{{\rm d}L} \over {{\rm d}t}} = K_{\rm a} \mathop \sum \limits_{\rm \sigma} b_{\rm a} ({\rm \sigma} ){\rm \Lambda} _{\rm \sigma} ^{\rm a} (x,e),$$

(15)$$\displaystyle{{{\rm d}X} \over {{\rm d}t}} = K_{\rm a} \mathop \sum \limits_{\rm \sigma} b_{\rm a} ({\rm \sigma} ){\rm \Upsilon} _{\rm \sigma} ^{\rm a} (x,e),$$

where

(16)$$K_{\rm a} = - \displaystyle{{3M_* R^3} \over {5{\rm \bar \rho} a^3}}. $$

The b _a(σ) factor is now:

(17)$$\eqalign{b_{\rm a} ({\rm \sigma} )& = \tilde p_2 ({\rm \sigma} )\sin {\rm \;} \left( {{\rm \sigma} {\rm \Delta} t_{\rm a} ({\rm \sigma} )} \right) = \tilde p_2 ({\rm \sigma} )\sin {\rm \;} 2{\rm \delta} _{\rm a} ({\rm \sigma} ) \cr & = |\tilde p_2 |\sin 2({\rm \delta} _{\rm a} ({\rm \sigma} ) + {\rm \pi} /2) = - |\tilde p_2 |\sin 2{\rm \delta} _{\rm a} ({\rm \sigma} ),} $$

where Δt _a is the atmosphere's delayed response to the stellar heat excitation, and δ_a is the corresponding geometric lag. The amplitude of the bulge, $\tilde p_2 $, is the second-order surface pressure variations (Chapman & Lindzen Reference Chapman and Lindzen1970):

(18)$$\tilde p_2 ({\rm \sigma} ) = {\rm i}\displaystyle{\gamma \over {\rm \sigma}} \tilde p_0 \left( {\nabla \cdot {\bf \upsilon} _{\rm \sigma} - \displaystyle{{{\rm \gamma} - 1} \over {\rm \gamma}} \displaystyle{{J_{\rm \sigma}} \over {gH_0}}} \right) = {\rm i}\displaystyle{{{\rm {\cal P}}_{\rm \sigma}} \over {\rm \sigma}}, $$

where γ=7/5 for a perfect diatomic gas, $\tilde p_0 $ is the mean surface pressure, υ_σ is the tidal winds velocity, J _σ is the amount of heat absorbed or emitted by unit of mass of air per unit time, and H ₀ is the scale height at the surface. The imaginary number in equation (18) causes the pressure variations to lead the star (i=e^iπ/2).

The coefficients Λ_σ^a and ${\rm \Upsilon} _{\rm \sigma} ^{\rm a} $ are also polynomials in the eccentricity, but different from their analogues for gravitational tides (equations (10) and (11)). However, for zero eccentricity they become equal, i.e., Λ_σ^a(e=0)=Λ_σ^g(e=0), and ${\rm \Upsilon} _{\rm \sigma} ^{\rm a} (e = 0) = {\rm \Upsilon} _{\rm \sigma} ^{\rm g} (e = 0)$. Once more, for a planet with moderate eccentricity, we can neglect terms in e ⁴ and greater, and obtain the following expansion for equations (14) and (15):

(19)$$\eqalign{\displaystyle{1 \over {K_{\rm a}}} \displaystyle{{{\rm d}L} \over {{\rm d}t}} =\ & b_{\rm a} (2{\rm \omega} + 3n)e^2 \displaystyle{{27} \over {32}}(1 - x)^4 \cr & + b_{\rm a} (2{\rm \omega} + 2n)\displaystyle{3 \over {32}}(1 - 3e^2 )(1 - x)^4 \cr & - b_{\rm a} (2{\rm \omega} + n)e^2 \displaystyle{3 \over {32}}(1 - x)^4 \cr & + b_{\rm a} (2{\rm \omega} )\displaystyle{3 \over 8}(1 + 5e^2 )(1 - x^2 )^2 \cr & - b_{\rm a} (2{\rm \omega} - n)e^2 \displaystyle{3 \over {32}}(1 + x)^4 \cr & + b_{\rm a} (2{\rm \omega} - 2n)\displaystyle{3 \over {32}}(1 - 3e^2 )(1 + x)^4 \cr & + b_{\rm a} (2{\rm \omega} - 3n)e^2 \displaystyle{{27} \over {32}}(1 + x)^4 \cr & + b_{\rm a} ({\rm \omega} + 3n)e^2 \displaystyle{{27} \over {16}}(1 - x^2 )(1 - x)^2 \cr & + b_{\rm a} ({\rm \omega} + 2n)\displaystyle{3 \over {16}}(1 - 3e^2 )(1 - x^2 )(1 - x)^2 \cr & - b_{\rm a} ({\rm \omega} + n)e^2 \displaystyle{3 \over {16}}(1 - x^2 )(1 - x)^2 \cr & + b_{\rm a} ({\rm \omega} )\displaystyle{3 \over 4}(1 + 5e^2 )x^2 (1 - x^2 ) \cr & - b_{\rm a} ({\rm \omega} - n)e^2 \displaystyle{3 \over {16}}(1 - x^2 )(1 + x)^2 \cr & + b_{\rm a} ({\rm \omega} - 2n) \displaystyle{3 \over {16}}(1 - 3e^2 )(1 - x^2 )(1 + x)^2 \cr & + b_{\rm a} ({\rm \omega} - 3n)e^2 \displaystyle{{27} \over {16}}(1 - x^2 )(1 + x)^2,} $$

and

(20)$$\eqalign{\displaystyle{1 \over {K_{\rm a}}} \displaystyle{{{\rm d}X} \over {{\rm d}t}} = & - b_{\rm a} (2{\rm \omega} + 3n)e^2 \displaystyle{{27} \over {32}}(1 - x)^4 \cr & - b_{\rm a} (2{\rm \omega} + 2n)\displaystyle{3 \over {32}}(1 - 3e^2 )(1 - x)^4 \cr & + b_{\rm a} (2{\rm \omega} + n)e^2 \displaystyle{3 \over {32}}(1 - x)^4 \cr & - b_{\rm a} (2{\rm \omega} - n)e^2 \displaystyle{3 \over {32}}(1 + x)^4 \cr & + b_{\rm a} (2{\rm \omega} - 2n)\displaystyle{3 \over {32}}(1 - 3e^2 )(1 + x)^4 \cr & + b_{\rm a} (2{\rm \omega} - 3n)e^2 \displaystyle{{27} \over {32}}(1 + x)^4 \cr & - b_{\rm a} ({\rm \omega} + 3n)e^2 \displaystyle{{27} \over 8}(1 - x^2 )(1 - x)^2 \cr & - b_{\rm a} ({\rm \omega} + 2n)\displaystyle{3 \over 8}(1 - 3e^2 )(1 - x^2 )(1 - x)^2 \cr & + b_{\rm a} ({\rm \omega} + n)e^2 \displaystyle{3 \over 8}(1 - x^2 )(1 - x)^2 \cr & - b_{\rm a} ({\rm \omega} - n)e^2 \displaystyle{3 \over 8}(1 - x^2 )(1 + x)^2 \cr & + b_{\rm a} ({\rm \omega} - 2n)\displaystyle{3 \over 8}(1 - 3e^2 )(1 - x^2 )(1 + x)^2 \cr & + b_{\rm a} ({\rm \omega} - 3n)e^2 \displaystyle{{27} \over 8}(1 - x^2 )(1 + x)^2 \cr & - b_{\rm a} (3n)e^2 \displaystyle{{81} \over {16}}(1 - x^2 )^2 \cr & - b_{\rm a} (2n)\displaystyle{9 \over {16}}(1 - 3e^2 )(1 - x^2 )^2 \cr & + b_{\rm a} (n)e^2 \displaystyle{9 \over {16}}(1 - x^2 )^2.} $$

Neglecting the tidal winds υ_σ, from expression (18) we have

(21)$${\rm {\cal P}}_{\rm \sigma} \propto J_{\rm \sigma} \propto r^{ - 2}. $$

As a consequence, in the computation of the averaged expressions for Λ_σ^a and ${\rm \Upsilon} _{\rm \sigma} ^{\rm a} $ we used a factor r^{− 5}, corresponding to the r^{− 2} from ${\rm {\cal P}}_{\rm \sigma} $ together with the contribution in r^{− 3} from the tidal potential (equation (13)).

Tidal models

The dissipation of the mechanical energy of tides in the planet's interior is responsible for the phase lags δ(σ). A commonly used dimensionless measure of tidal damping is the quality factor Q, defined as the inverse of the ‘specific’ dissipation and related to the phase lags by Efroimsky (Reference Efroimsky2012) (equation (141))

(22)$$Q^{ - 1} ({\rm \sigma} ) = \displaystyle{{\Delta E} \over {2{\rm \pi} E}} = |\sin 2{\rm \delta} _{\rm g} ({\rm \sigma} )|,$$

where E is the total tidal energy stored in the planet, and ΔE is the energy dissipated per cycle. We can rewrite expression (12) as:

(23)$$b_{\rm g} ({\rm \sigma} ) = {\rm sign}({\rm \sigma} )\displaystyle{{k_2} \over {Q({\rm \sigma} )}}.$$

The present Q value for the planets in the Solar System can be estimated from orbital measurements, but as rheology of the planets is badly known, the exact dependence of b _g(σ) on the tidal frequency σ is unknown. Many different authors have studied the problem and several models for b _g(σ) have been developed so far, from the simplest ones to the more complex (for a review see Efroimsky & Williams Reference Efroimsky and Williams2009). The huge problem in validating one model better than the others is the difficulty to compare the theoretical results with the observations, as the effect of tides are very small and can only be detected efficiently after long periods of time. Therefore, here we only describe the most commonly models that are used (Fig. 2).

Fig. 2. Commonly used models for frequency dependence of tides: visco-elasctic (equation (24), red), viscous or linear (equation (26), blue), constant – Q (equation (27), green), and an interpolated model between the linear and the constant one (equation (75), dashed line).

The visco-elastic model: In this model, Darwin (Reference Darwin1908) assumed that the planet behaves like a Maxwell solidFootnote ⁷ of constant density ρ, and found:

(24)$$b_{\rm g} ({\rm \sigma} ) = k_{\rm f} \displaystyle{{{\rm \tau} _{\rm b} - {\rm \tau} _{\rm a}} \over {1 + ({\rm \tau} _{\rm b} {\rm \sigma} )^2}} {\rm \sigma}, $$

where τ_a=υ_e/μ_e and τ_b are the time constants for damping of the body tides,

(25)$${\rm \tau} _{\rm b} = {\rm \tau} _{\rm a} (1 + 19{\rm \mu} _{\rm e} R/2Gm{\rm \rho} ).$$

The Maxwell solid model is usually accepted as a good approximation of the planet's response to tidal perturbations, although more complete visco-elastic models exist (for a review see Henning et al. Reference Henning, O'Connell and Sasselov2009). However, it depends on many uncertain parameters and it is not of practical use when we aim to do simple approximations and considerations on the tidal evolution. In addition, while it can be appropriate to simultaneously describe the elastic and viscous response of the solid body, it is questionable whether it is still valid for the atmospheric deformation is response to thermal gradients.

The viscous or linear model: in this model, it is assumed that the response time delay to the perturbation is independent of the tidal frequency, i.e., the position of the ‘maximal tide’ is shifted from the sub-stellar point by a constant time-lag Δt _g (e.g, Mignard Reference Mignard1979). As usual we have σΔt _g≪1, the viscous model becomes linear:

(26)$$b_{\rm g} ({\rm \sigma} ) = k_2 \sin \,({\rm \sigma} {\rm \Delta} t_{\rm g} ) \approx k_2 \,{\rm \sigma} {\rm \Delta} t_{\rm g}. $$

The viscous model is a particular case of the visco-elastic model and is specially adapted to describe the behaviour of planets in slow rotating regimes (ω~n).

The constant-Q model: for the Earth, Q changes by slightly more than an order of magnitude between the Chandler wobble period (about 440 days) and seismic periods of a few seconds (e.g., Anderson & Minster Reference Anderson and Minster1979; Karato & Spetzler Reference Karato and Spetzler1990). Thus, it is also common to treat the specific dissipation as independent of frequency:

(27)$$b_{\rm g} ({\rm \sigma} ) \approx {\rm sign}({\rm \sigma} )k_2 /Q.$$

The constant-Q model can be used for long periods of time where the tidal dissipation does not change much, as is the case for moderately rotating planets. However, for slow rotating planets, the constant-Q model is not appropriate as it gives rise to discontinuities for σ=0.

Core–mantle friction

We assume that the Earth-sized planets have internal structures similar to that of our planet, that is, these planets are composed of a solid mantle and a liquid core. We suppose that the liquid core is inviscid, incompressible and homogenous. We also suppose that the internal structure of the planets is unchanged in time, since the core formation time (or condensation velocity) is poorly known.

If there is slippage between the liquid core and the mantle, an additional source of dissipation of rotational energy results from friction occurring at the core–mantle boundary. Indeed, because of their different shapes and densities, the core and the mantle do not have the same dynamical ellipticity and the two parts tend to precess at different rates (Poincaré Reference Poincaré1910). This tendency is more or less counteracted by different interactions produced at their interface: the torque N of non-radial inertial pressure forces of the mantle over the core provoked by the non-spherical shape of the interface; the torque of the viscous (or turbulent) friction between the core and the mantle; and the torque of the electromagnetic friction, caused by the interaction between electrical currents of the core and the bottom of the magnetized mantle.

Poincaré's (Reference Poincaré1910) classic study shows that the core responds to the precession with a rotation velocity ω_c, whose vector is inclined by a small angle χ with respect to the mantle velocity vector ω=ωk (Rochester Reference Rochester1976; Correia Reference Correia2006):

(28)$${\bf \delta} = {\bf \omega} - {\rm \omega} _{\rm c}, $$

and

(29)$$\sin {\rm \chi} = \left| {\displaystyle{{{\bf \delta} \times {\bf k}} \over {{\rm \omega} _{\rm c}}}} \right|.$$

Sasao et al. (Reference Sasao, Okubo and Saito1980) demonstrated that the pressure inertial torque may be expressed in a general way by:

(30)$${\bi N} = {\rm \omega} _{\rm c} \times {\bi L}_{\rm c} - {\bi P}_{\rm c}, $$

where P_c is the precession torque over the core and L_c its angular moment:

(31)$${\bi L}_{\rm c} = {{\bf \tilde { {\hskip2\cal I}}}}_{\rm c} \cdot {\bf \omega} _{\rm c}, $$

with

(32)$${{ \bf\tilde {{\hskip3 \cal I}}}}_{\rm c} = \left( {\matrix{ {A_{\rm c}} & 0 & 0 \cr 0 & {A_{\rm c}} & 0 \cr 0 & 0 & {C_{\rm c}} \cr}} \right).$$

C _c and A _c are the principal inertial moments of the core. The core's dynamical ellipticity is then given by E _c=(C _c−A _c)/C _c. The lower part of the mantle has irregularities that prevent the border between the two layers from being a perfect ellipsoid. According to Hide (Reference Hide1969) the height of these bumps may reach a few kilometres. Thus, for the dynamical ellipticity of the mantle, we establish that E _c=E _ch+δE _c, that is, E _c has a hydrostatic component and a non-hydrostatic component. For Earth, Herring et al. (Reference Herring, Gwinn and Shapiro1986) provide the value δE _c=1.2×10⁻⁴.

The two types of friction torques (viscous and electromagnetic) depend on the differential rotation between the core and the mantle and can be expressed by a single effective friction torque, Φ. As a general expression for this torque, we adopt (Rochester Reference Rochester1976; Mathews & Guo Reference Mathews and Guo2005)

(33)$${\bf \Phi} = - ({\rm \kappa} + {\rm \kappa '}\,{\bf k} \times )\,{\bf \delta}, $$

where κ and κ′ are effective coupling parameters. Both parameters account for viscous and electromagnetic stresses at the core–mantle interface that can be written as: κ=κ_vis+κ_em and κ′=κ′_vis+κ′_em. Estimations for these coefficients can be found in the works of Mathews & Guo (Reference Mathews and Guo2005) and Deleplace & Cardin (Reference Deleplace and Cardin2006). In the simplified case of no magnetic field, the coupling parameters are only given by the viscous friction contributions, which can be simplified as (Noir et al. Reference Noir, Cardin, Jault and Masson2003; Mathews & Guo Reference Mathews and Guo2005):

(34)$${\rm \kappa} _{{\rm vis}} = 2.62\sqrt {{\rm \nu} |{\rm \omega} |} /R_{\rm c} \quad {\rm and}\quad {\rm \kappa '}_{{\rm vis}} = 0.259\sqrt {{\rm \nu} |{\rm \omega} |} /R_{\rm c}, $$

where R _c is the core radius and ν the kinematic viscosity, which is poorly known. Even in the case of the Earth, the uncertainty in ν covers about 13 orders of magnitude (Lumb & Aldridge Reference Lumb and Aldridge1991), the best estimate so far being ν≃10⁻⁶ m² s⁻¹ (Gans Reference Gans1972; Poirier Reference Poirier1988; Wijs et al. Reference Wijs, Kresse and Vočadlo1998).

The contributions for spin variations of a planet may be obtained writing the equations of angular momentum conservation for the core and mantle:

(35)$$\displaystyle{{{\rm d}{\bi L}_{\rm m}} \over {{\rm d}t}} = {\bi P}_{\rm m} - {\bi N} - {\bf \Phi}, $$

(36)$$\displaystyle{{{\rm d}{\bi L}_{\rm c}} \over {{\rm d}t}} = {\bi P}_{\rm c} + {\bi N} + {\bf \Phi}, $$

where the P_m is the precession torque over the mantle and L_m is its angular momentum.

A general formulation for the equations of motion, valid for the both fast and slow rotation regimes of the planet, can be written as (Correia Reference Correia2006):

(37)$$\displaystyle{{{\rm d}L} \over {{\rm d}t}} \simeq - \displaystyle{{{\rm \kappa} A_{\rm c} C_{\rm c} {\rm \omega \alpha} ^2 \cos ^2 {\rm \varepsilon} \sin ^2 {\rm \varepsilon}} \over {(C_{\rm c} E_{\rm c} {\rm \omega} )^2 + {\rm \kappa} ^2}}, $$

(38)$$\hskip-98pt \displaystyle{{{\rm d}X} \over {{\rm d}t}} \simeq 0.$$

Obliquity evolution

Until now, we have been expressing the variations of the spin in Andoyer's variables. Despite their practical use, these variables do not give a clear view of the obliquity variation. Since x=cos ε=X/L, one obtains:

(39)$$\displaystyle{{{\rm d}x} \over {{\rm d}t}} = \displaystyle{1 \over L}\left( {\displaystyle{{{\rm d}X} \over {{\rm d}t}} - x\displaystyle{{{\rm d}L} \over {{\rm d}t}}} \right).$$

For tidal effects (τ=g or τ=a), we express dx/dt using the eccentricity series for dL/dt and dX/dt (equations (10, 11) and (19, 20)):

(40)$$\eqalign{\displaystyle{{{\rm d}x} \over {{\rm d}t}} = & - \displaystyle{{K_{\rm \tau}} \over {\rm \omega}} \mathop \sum \limits_{\rm \sigma} b_{\rm \tau} ({\rm \sigma} )({\rm \Upsilon} _{\rm \sigma} - x{\rm \Lambda} _{\rm \sigma} ) \cr = & K_{\rm \tau} \displaystyle{{1 - x^2} \over {\rm \omega}} \mathop \sum \limits_{\rm \sigma} b_{\rm \tau} ({\rm \sigma} ){\rm \Theta} _{\rm \sigma} (x,e).} $$

We thus have

(41)$$\displaystyle{{{\rm d}x} \over {{\rm d}t}} \propto \displaystyle{{1 - x^2} \over {\rm \omega}} \displaystyle{{{\rm d}{\rm \omega}} \over {{\rm d}t}},$$

meaning that the obliquity variations are smaller than the rotation rate variations for initially fast rotating planets, and that ${\rm \dot \varepsilon} \propto - \sin {\rm \varepsilon} $.

For the core–mantle friction effect, the variation of ε is easily computed from expression (39), since dX/dt≃0 (equation (38)).

(42)$$\displaystyle{{{\rm d}x} \over {{\rm d}t}} \simeq - \displaystyle{x \over L}\displaystyle{{{\rm d}L} \over {{\rm d}t}} = - \displaystyle{x \over {\rm \omega}} \displaystyle{{{\rm d}{\rm \omega}} \over {{\rm d}t}}$$

It also follows that

(43)$$X = cte.\quad \Rightarrow \quad x\,{\rm \omega} = cte.$$

Gravitational tides alone

From expression (18), we notice that for the initial stages of the evolution (σ≫n) atmospheric tides are weak. The same is true for core–mantle friction (equation (37)). As a consequence, gravitational tides dominate the spin evolution for fast rotating rates, thermal tides and core–mantle friction only playing a role for slow rotations (Correia et al. Reference Correia, Laskar and Néron de Surgy2003). Therefore, most studies on the spin evolution only consider gravitational tides. For a better comparison with our study, we recall here the main consequences of this effect alone.

If we assume a fast initial rotation, planets spend most of their evolution in the ω≫n regime. In this case, dω/dt evolves independently of the dissipative model, because all the terms in expression (10) have the same sign. Thus, in this regime gravitational tides always decrease the rotation rate, since K _g<0 (equation (9)).

When the planet arrives in the slow rotation regime (ω~n), some of the terms in equation (10) become negative and we can no longer generalize our conclusions to all dissipative models. However, in the vicinity of tidal frequencies near zero (σ≈0), except for the constant – Q model, all dissipative models can be made linear (Fig. 2). Then, adopting the linear dissipation model described in the ‘Tidal models’ subsection, we can write

(44)$$b_{\rm g} ({\rm \sigma} ) \simeq k_2 {\rm \sigma} {\rm \Delta} t_{\rm g} = \displaystyle{{k_2} \over {Q_n}} \displaystyle{{\rm \sigma} \over {2n}},$$

where Q _n⁻¹=2nΔt _g is the specific dissipation factor for small frequencies. For small eccentricity, the highest tidal frequency is σ=2ω+3n (equation (19)), so this approximation is validFootnote ⁸ for Q _n≳5.

In the limit of slow rotation rates, we can then simplify expression (10):

(45)$$\displaystyle{{{\rm d}{\rm \omega}} \over {{\rm d}t}} = - K_0 \left[ {{\rm \omega} \left( {\displaystyle{{1 + x^2} \over 2}} \right) - n(1 + 6e^2 )x} \right]\;, $$

where

(46)$$K_0 = - \displaystyle{{3K_{\rm g} k_2 \Delta t_{\rm g}} \over C}\left( {1 + \displaystyle{{15} \over 2}e^2} \right)\;. $$

Thus, for each obliquity there is an equilibrium value for the rotation rate ω_e, obtained when dω/dt=0:

(47)$$\displaystyle{{{\rm \omega} _{\rm e}} \over n} = (1 + 6e^2 )\displaystyle{{2x} \over {1 + x^2}}. $$

For ω>ω_e the rotation rate decreases, whereas for ω<ω_e it increases.

Equilibrium damping time

In order to check if the spin of a given planet is still evolving or if it can be found at its equilibrium rotation, we can compute an evolutionary characteristic time and then compare it to the age of the host star. Since gravitational tides are the dominant effect for the most part of the evolution, we can use this effect alone to estimate the characteristic time τ_eq needed to reach the equilibrium rotation. Thus, from expression (45) we get:

(48)$${\rm \tau} _{{\rm eq}} \sim \displaystyle{1 \over {K_0}} = \displaystyle{{GC} \over {3k_2 \Delta t_{\rm g} R^5}} \displaystyle{1 \over {n^4}} = \displaystyle{{ma^6} \over {9GM_*^2 k_2 \Delta t_{\rm g} R^3}}, $$

with C≈m R ²/3 and n ⁴=G ²M _*²/a ⁶. In Table 1, we computed the characteristic times for existing Earth-sized exoplanets. We notice that, except for GJ 667C g, all known planets have a τ_eq lower than the age of the system, and so it is believed that they have already reached an equilibrium rotation state.

Table 1. Characteristics and equilibrium rotation rates of Earth-sized planets with masses lower than 10 M _⊕ (see text for notations)

^a Using Equation (48) with k ₂=1/3 and Δt _g=640 s (Earth's values). References: 1) Rappaport et al. (Reference Rappaport, Sanchis-Ojeda, Rogers, Levine and Winn2013b); 2) Charpinet et al. (Reference Charpinet2011); 3) Rappaport et al. (Reference Rappaport, Barclay and DeVore2013a); 4) Muirhead et al. (Reference Muirhead, Johnson and Apps2012); 5) Callegari & Rodríguez (Reference Callegari and Rodríguez2013); 6) Dumusque et al. (Reference Dumusque, Pepe and Lovis2012); 7) Xie (Reference Xie2013); 8) Hadden & Lithwick (Reference Hadden and Lithwick2013); 9) Chaplin et al. (Reference Chaplin, Sanchis-Ojeda and Campante2013); 10) Sanchis-Ojeda et al. (Reference Sanchis-Ojeda, Rappaport and Winn2013); 11) Pepe et al. (Reference Pepe, Cameron and Latham2013); 12) Lopez & Fortney (Reference Lopez and Fortney2013); 13) Lissauer et al. (Reference Lissauer, Fabrycky and Ford2011); 14) Forveille et al. (Reference Forveille, Bonfils and Delfosse2011); 15) Pepe et al. (Reference Pepe, Lovis and Ségransan2011); 16) Anglada-Escudé et al. (Reference Anglada-Escudé, Tuomi and Gerlach2013); 17) Fressin et al. (Reference Fressin, Torres and Rowe2012); 18) Tuomi et al. (Reference Tuomi, Anglada-Escudé and Gerlach2013); 19) Borucki et al. (Reference Borucki, Agol and Fressin2013); 20) Jontof-Hutter et al. (Reference Jontof-Hutter, Lissauer, Rowe and Fabrycky2013); 21) Howard et al. (Reference Howard, Johnson and Marcy2011); 22) Carter et al. (Reference Carter, Agol and Chaplin2012); 23) Anglada-Escudé & Tuomi (Reference Anglada-Escudé and Tuomi2012); 24) Batalha et al. (Reference Batalha, Borucki and Bryson2011); 25) Gilliland et al. (Reference Gilliland2013); 26) Vogt et al. (Reference Vogt, Wittenmyer and Butler2010) 27) Steffen et al. (Reference Steffen, Fabrycky and Agol2013); 28) Delfosse et al. (Reference Delfosse, Bonfils and Forveille2013); 29) Rivera et al. (Reference Rivera, Butler and Vogt2010a); 30) Harpsøe et al. (Reference Harpsøe, Hardis and Hinse2013); 31) Ofir et al. (Reference Ofir, Dreizler, Zechmeister and Husser2013); 32) Lo Curto et al. (Reference Lo Curto, Mayor and Benz2010); 33) Rivera et al. (Reference Rivera, Laughlin and Butler2010b); 34) Cochran et al. (Reference Cochran2011); 35) Bonfils et al. (Reference Bonfils, Gillon and Forveille2011); 36) Holman et al. (Reference Holman, Fabrycky and Ragozzine2010); 37) Torres et al. (Reference Torres, Fressin and Batalha2011); 38) Léger et al. (Reference Léger, Rouan and Schneider2009); 39) Bouchy et al. (Reference Bouchy, Mayor and Lovis2009); 40) Henry et al. (Reference Henry, Howard, Marcy, Fischer and Johnson2011); 41) Endl et al. (Reference Endl, Robertson and Cochran2012); 42) Queloz et al. (Reference Queloz, Bouchy and Moutou2009); 43) Forveille et al. (Reference Forveille, Bonfils and Delfosse2009); 44) Lo Curto et al. (Reference Lo Curto, Mayor and Benz2013); 45) Gautier et al. (Reference Gautier, Charbonneau and Rowe2012); 46) Jenkins et al. (Reference Jenkins, Tuomi, Brasser, Ivanyuk and Murgas2013); 47) Howard et al. (Reference Howard, Johnson and Marcy2009); 48) Lovis et al. (Reference Lovis, Mayor and Pepe2006); 49) Howell et al. (Reference Howell, Rowe and Bryson2012); 50) Weiss & Marcy (Reference Weiss and Marcy2013); 51) Pepe et al. (Reference Pepe, Correia and Mayor2007).

Thermal tides effect

We now include in our analysis the effect from thermal atmospheric tides to the spin evolution (subsection ‘Thermal atmospheric tides’). As shown earlier, the averaged variation of the spin can be expressed by equation (10) for gravitational tides and by equation (19) for atmospheric tides. Considering that the planet's obliquity is small (ε~0), we can neglect terms of order ε² or higher in equations (10) and (19). In this case, the variation of the spin caused by gravitational tides is then (Correia et al. Reference Correia, Levrard and Laskar2008):

(49)$$\eqalign{ T_{\rm g} =\ &\left. {\displaystyle{{{\rm d}L} \over {{\rm d}t}}} \right]_{\rm g} = \displaystyle{3 \over 2}K_{\rm g} \left[ {(1 - 5e^2 )b_{\rm g} (2{\rm \omega} - 2n)} \right. \cr & \left. { + \displaystyle{1 \over 4}e^2 b_{\rm g} (2{\rm \omega} - n) + \displaystyle{{49} \over 4}e^2 b_{\rm g}\right. \left. (2{\rm \omega} - 3n)} \big].} $$

In the same way, for thermal atmospheric tides we obtain

(50)$$\eqalign{ T_{\rm a} =\ &\left. {\displaystyle{{{\rm d}L} \over {{\rm d}t}}} \right]_{\rm a} = \displaystyle{3 \over 2}K_{\rm a} \left[ {(1 - 3e^2 )b_{\rm a} (2{\rm \omega} - 2n)} \right. \cr & \left. { - e^2 b_{\rm a} (2{\rm \omega} - n) + 9e^2 b_{\rm a} (2{\rm \omega} - 3n)} \right]\;.} $$

As the time-lag from gravitational and atmospheric tides is poorly known, in this section we use again the viscous model. Since σΔt is usually small, then:

(51)$$b_{\rm g} ({\rm \sigma} ) \simeq k_2 {\rm \sigma} {\rm \Delta} t_{\rm g} \quad {\rm and}\quad b_{\rm a} ({\rm \sigma} ) \simeq - |\tilde p_2 |{\rm \sigma} {\rm \Delta} t_{\rm a}. $$

For atmospheric tides it is also necessary to consider the response of the surface pressure variations to tidal frequency (equation (18)). We use the ‘heating at the ground’ model, described by Dobrovolskis & Ingersoll (Reference Dobrovolskis and Ingersoll1980). It is supposed that all the stellar flux absorbed by the ground, F _s, is immediately deposited in a thin layer of atmosphere at the surface. The heating distribution is then written as a delta-function just above the ground ($J_{\rm \sigma} = gF_{\rm s} /\tilde p_0 $). Besides the good agreement with the observations, this approximation is justified because tides in the upper atmosphere are decoupled from the ground by the disparity between their rotation rates. Neglecting υ_σ over the thin heated layer, equation (18) becomes:

(52)$${\rm {\cal P}}_{\rm \sigma} = F_{\rm s} /(8H_0 ) \propto L_* /a^2, $$

where L _* is the star luminosity and a the semi-major axis.

The average evolution of the rotation rate can then be obtained adding the effects of both tidal torques acting on the planet:

(53)$$\displaystyle{{{\rm d}L} \over {{\rm d}t}} = C{\rm \dot \omega} = (T_{\rm g} + T_{\rm a} ).$$

Substituting in previous expression T _g and T _a given by equations (49) and (50) together with the viscous dissipative model (equation (51)), we haveFootnote ⁹:

(54)$$\eqalign{ \displaystyle{{{\rm \dot \omega}} \over { - K_0}} =\ & {\rm \omega} - (1 + 6e^2 )n - {\rm \omega} _{\rm s} \left[ {\left( {1 - \displaystyle{{21} \over 2}e^2} \right){\rm sign}({\rm \omega} - n)} \right. \cr & \left. { - e^2 {\rm sign}(2{\rm \omega} - n) + 9e^2 {\rm sign}(2{\rm \omega} - 3n)} \bigg]} $$

with

(55)$${\rm \omega} _{\rm s} = \displaystyle{{K_{\rm a} F_{\rm s} \Delta t_{\rm a}} \over {16H_0 K_{\rm g} k_2 \Delta t_{\rm g}}} \propto \displaystyle{{L_* R} \over {M_* m}}a.$$

Equilibrium final states for the rotation rate

An equilibrium final state occurs when ${\rm \dot \omega} = 0$. From equation (53) this happens when T _g=−T _a. When e=0 (see Correia & Laskar Reference Correia and Laskar2001),

(56)$$f\,({\rm \omega} - n) = - \displaystyle{{T_{\rm g}} \over {T_{\rm a}}} = - \displaystyle{{K_{\rm g} b_{\rm g} (2{\rm \omega} - 2n)} \over {K_{\rm a} b_{\rm a} (2{\rm \omega} - 2n)}} = 1\;, $$

where f(x) is an even function of x. For the dissipative model in use, we can assume that f(x) is monotonic close to equilibrium, and so we haveFootnote ¹⁰:

(57)$$|{\rm \omega} - n| = f^{ - 1} (1) = {\rm \omega} _{\rm s}. $$

This means that there are two final possibilities for the equilibrium rotation of the planet:

(58)$${\rm \omega} ^ \pm = n \pm {\rm \omega} _{\rm s}. $$

If ω_s<n, all final states correspond to prograde final rotation rates (Fig. 3(a) and (b)). Retrograde rotation appears only if the planet evolves to the ω⁻ state together with ω_s>n. As we can see in Table 1, this is the case of planet Venus (Fig. 3(c)).

Fig. 3. Variation of ${\rm \dot \omega} $ with ω/n (equation (54)) for (a) ω_s/n=0.05, (b) ω_s/n=0.55 and (c) ω_s/n=1.92, using different eccentricities (e=0.0, 0.1, 0.2). The equilibrium rotation rates are given by ${\rm \dot \omega} $=0 and the arrows indicate whether it is a stable or unstable equilibrium position.

For moderate values of the eccentricity and using equations (49) and (50), equation (56) can be rewritten as:

(59)$$f\,({\rm \omega} - n) = 1 - e^2 \left( {2 + \displaystyle{{g({\rm \omega} )} \over {b_{\rm a} (2{\rm \omega} - 2n)}}} \right),$$

where

(60)$$\eqalign{g({\rm \omega} ) = & \displaystyle{{K_{\rm g}} \over {4K_{\rm a}}} [b_{\rm g} (2{\rm \omega} - n) + 49b_{\rm g} (2{\rm \omega} - 3n)] \cr & - b_{\rm a} (2{\rm \omega} - n) + 9b_{\rm a} (2{\rm \omega} - 3n).} $$

When we compute the inverse function of previous equation and use ω_s from expression (57) we get:

(61)$$|{\rm \omega} - n| = f^{ - 1} (1) - e^2 \left( {2 + \displaystyle{{g({\rm \omega} )} \over {b_{\rm a} (2{\rm \omega} - 2n)}}} \right)\left. {\displaystyle{{\partial f^{ - 1}} \over {\partial x}}} \right|_{x = 1}, $$

which means that we now have four final possibilities for the equilibrium rotation of the planet:

(62)$${\rm \omega} _{1,2}^ \pm = n \pm {\rm \omega} _{\rm s} + e^2 {\rm \delta} _{1,2}^ \pm, $$

with

(63)$${\rm \delta} _{1,2}^ \pm = \left( {2 + \displaystyle{{g({\rm \omega} )} \over {|b_{\rm a} (2{\rm \omega} - 2n)|}}} \right)\left. {\displaystyle{{\partial f^{ - 1}} \over {\partial x}}} \right|_{x = 1} \;. $$

Adopting the dissipation model from equation (51) and using expression (54), we find that ω₂⁺ occurs when ω>3n/2 and

(64)$$\eqalign{{\rm \omega} _2^ + =\ & (1 + 6e^2 )n + {\rm \omega} _{\rm s} \left[ {1 + e^2 \left( {9 - \displaystyle{{21} \over 2} - 1} \right)} \right] \cr =\ & n + {\rm \omega} _{\rm s} + e^2 \left( {6n - \displaystyle{5 \over 2}{\rm \omega} _{\rm s}} \right).} $$

Similarly, we find the ω₁⁺ state when n<ω<3n/2

(65)$${\rm \omega} _1^ + = n + {\rm \omega} _{\rm s} + e^2 \left( {6n - \displaystyle{{41} \over 2}{\rm \omega} _{\rm s}} \right),$$

when n/2<ω<n we have the ω₁⁻ state

(66)$${\rm \omega} _1^ - = n - {\rm \omega} _{\rm s} + e^2 \left( {6n + \displaystyle{1 \over 2}{\rm \omega} _{\rm s}} \right),$$

and finally for ω<n/2 we have

(67)$${\rm \omega} _2^ - = n - {\rm \omega} _{\rm s} + e^2 \left( {6n + \displaystyle{5 \over 2}{\rm \omega} _{\rm s}} \right).$$

Then, expression (63) can be rewritten as:

(68)$${\rm \delta} _{1,2}^ - = 6n + \left( {\displaystyle{3 \over 2} \pm 1} \right){\rm \omega} _{\rm s} $$

and

(69)$${\rm \delta} _{1,2}^ + = 6n - \left( {\displaystyle{{23} \over 2} \pm 9} \right){\rm \omega} _{\rm s}. $$

As the set of final states ω_1,2^± must also verify

(70)$${\rm \omega} _2^ - \lt n/2 \lt {\rm \omega} _1^ - \lt n \lt {\rm \omega} _1^ + \lt 3n/2 \lt {\rm \omega} _2^ +, $$

in general, and depending on the values of ω_s and e, these four equilibrium rotation states cannot coexist all simultaneously. In particular, when we adopt the viscous tidal model (equation (51)), the final states ω₁⁻ and ω₁⁺ never coexist with ω₂⁻. At most three different equilibrium states are therefore possible, obtained when ω_s/n is close to 1/2, or more precisely, when

(71)$$\displaystyle{{1/2 - 6e^2} \over {1 - 5e^2 /2}} \lt \displaystyle{{{\rm \omega} _{\rm s}} \over n} \lt \displaystyle{{1/2 - 6e^2} \over {1 - 41e^2 /2}}.$$

Conversely, we have that ω₁⁺ is the single final state that exists whenever ω_s/n<6e ²(1+e ²/2).

In Fig. 3, we plot some examples for the rotation rate evolution for some ω_s and eccentricity values. We see that when ω_s/n is close to zero only the final equilibrium state ω₁⁺ is possible (Fig. 3(a)). For ω_s/n=1.92, the Venus value, two final equilibrium states are possible: the ω₂⁻ and the ω₂⁺ (Fig. 3(c)). However, when ω_s/n=0.55 we still have two possible states for e=0.0 (ω₁⁻ and ω₁⁺), but there are three possible final states for e>0.1 (ω₁⁻,ω₁⁺ and ω₂⁺) (Fig. 3(b)).

In Fig. 4, we plot the possible equilibrium rotation states as a function of ω_s/n for eccentricities of e=0.0, 0.1 and 0.2. Depending on ω_s/n and the eccentricity values, the planet may have one, two or three possibilities to evolve. If ω_s/n is close to 0.5 the planet has more possible final states (equation (71)): if the eccentricity is zero, the planet has two possible final states, but if the eccentricity of the planet is 0.1 or 0.2 there are three possibilities. We can also note that for e=0, retrogade states appear when ω_s/n⩾1, but for e=0.1 and e=0.2 the rotation rate must be slightly higher. More exactly, it is needed that ω_s/n⩾(1+17e ²/2), i.e., ω_s/n⩾1.09 for e=0.1, and ω_s/n⩾1.34 for e=0.2.

Fig. 4. Equilibrium positions of the rotation rate as a function of the ratio ω_s/n for three different values of the eccentricity: (a) e=0.0, (b) e=0.1 and (c) e=0.2. The solid red line corresponds to the ω₁⁺ state, the dotted red to the ω₁⁻ state, the solid blue line to the ω₂⁺ state, and the dotted blue line to the ω₂⁻ state.

Application to already known exoplanets

Although we have witnessed in the last few years an increase in the discovery of Earth-sized exoplanets, the physical data for these planets are still scarce. Hence, we can only do some assumptions that allow us to infer some constraints on their final spin evolution. We assume here that these planets are rocky with a dense atmosphere like Venus (Alibert et al. Reference Alibert, Baraffe and Benz2006; Rafikov Reference Rafikov2006). We can call these Earth-sized exoplanets as V-type planets. Using the empirical mass–luminosity relation L _*∝M _*⁴ (Cester et al. Reference Cester, Ferluga and Boehm1983) and the mass–radius relation for terrestrial planetsFootnote ¹¹R ∝ m ^0.274 (Sotin et al. Reference Sotin, Grasset and Mocquet2007), we get from equation (55):

(72)$$\displaystyle{{{\rm \omega} _{\rm s}} \over n} \propto \displaystyle{{L_* R} \over {M_* m}}\displaystyle{a \over n} = k(aM_* )^{2.5} m^{ - 0.726}. $$

k is a coefficient of proportionality containing the parameters that we cannot constrain for these planets (H ₀, F _s, k ₂, Δt _g and Δt _a). Assuming that all unknown properties of the Earth-sized planets are the same as for Venus (2π/ω_s=116.1 day, and so ω_s/n=1.9255), we computeFootnote ¹²:

(73)$$k = \displaystyle{{{\rm \omega} _{\rm s}} \over n}(aM_* )^{ - 2.5} m^{0.726} = 3.73({\rm AU}\;M_ \odot )^{ - 2.5} \,m_ \oplus ^{0.726} $$

Using the previous equation, we can estimate the ω_s/n ratio for all known Earth-sized planets (Table 1).

As shown in equation (72), ω_s/n is a function of the product aM _*. To see how the equilibrium rotation rate states depend on these parameters, in Fig. 5 we plot the number of values of the equilibrium rotation states as a function of aM _*, for eccentricities e=0.0, 0.1 and 0.2. We can compare Fig. 5 with Table 1, where the planet's parameters and the possible final equilibrium states are presented. The actual rotation state of the planet Venus corresponds to the ω₂⁻ equilibrium state. For the remaining exoplanets here studied, only one equilibrium rotation state is possible, the ω₁⁺ state. There is only one exception, planet HD 40307 f, which has two possible equilibrium states (ω₁⁺ and ω₁⁻), but their values are too close to be distinguished in a future direct observation of the planet.

Fig. 5. Equilibrium positions of the rotation rate as a function of the product aM _* for three different values of the eccentricity: (a) e=0.0, (b) e=0.1 and (c) e=0.2. The solid red line corresponds to the ω₁⁺ state, the dotted red line to the ω₁⁻ state, the solid blue line to the ω₂⁺ state, and the dotted blue line to the ω₂⁻ state.

The aM _* value for all planets listed in Table 1 is small as a consequence of the present limitations of the detection techniques: at present we are only able to detect Earth-sized planets orbiting very close to low-mass stars, for which gravitational tides dominate.

Generalization to Earth-sized planets

With the continuous increase in the precision of detection methods, we expect to soon find planets orbiting Sun-like stars as distant from them as the Earth or Venus are from the Sun (e.g., D'Odorico & CODEX/ESPRESSO Team Reference D'Odorico2007). Foreseeing this situation, we applied the previous calculations to long-period Earth-sized planets that may eventually be found in the HZ of these systems.

Assuming that the planet's mass and radius are the same of Venus and a host star with the mass of the Sun, we made calculations for three different values of the eccentricity (e=0.0, 0.1 and 0.2), varying for each one the value of the semi-major axis a. As expected, more possible equilibrium rotation states appear (Table 2). For the same eccentricity we note that the difference in the rotation period of the several equilibrium states is enlarged by the distance to the star. We can also see that for higher values of aM _* retrograde states may be possible. Setting ω₂⁻<0 (equation (67)) and substituting ω_s/n by equation (72), we may find the value of rotation rate needed for a retrograde state:

(74)$$aM_* \geqslant \left[ {\displaystyle{1 \over {3.73}}\left( {1 + \displaystyle{{17} \over 2}e^2} \right)} \right]^{0.4}. $$

Table 2. Characteristics and equilibrium rotation rates of Venus-like planets orbiting a Sun-like star (M_*=M_⊙)

Thus, for zero eccentricity a retrograde state is possible when aM _*≳0.59, for e=0.1 a possible retrograde states appears when aM _*≳0.61, whereas for an eccentricity of e=0.2, we need aM _*≳0.66 for retrograde rotation. This is also confirmed in Fig. 5, where we can see that planets with eccentricities of e=0.1 and 0.2 have a retrograde equilibrium state for an aM _* product higher than 0.6. Table 2 and Fig. 5 are also in agreement on the number of states when aM _* is around 0.45 AU M _⊙. Although only two states can be found when the eccentricity is zero, for eccentricities of 0.1 and 0.2 we find three equilibrium states.

An interpolated dissipative model

During the spin evolution of the planets, the tidal frequency varies and so does the dissipation factors b(σ) (equations (12) and (17)). Because the rheology of the planets is poorly known, the exact dependence of the dissipation on the tidal frequency is unknown. In the ‘Tidal models’ subsection, we described the most commonly used dissipation models (Fig. 2). However, the visco-elastic model seems inappropriate for atmospheric tides, whereas the viscous one is not realistic for σ≫n, and the constant – Q for σ=n. More complex models exist (e.g., Efroimsky & Williams Reference Efroimsky and Williams2009), but for simplicity we adopt here an interpolated model that behaves like the linear model for small values of σ and like the constant – Q model when σ increases (e.g., Ćuk & Stewart Reference Ćuk and Stewart2012). This model is also in agreement with the analytical simplifications from the ‘Dynamical analysis’ section. It is also very close to the models obtained by Remus et al. (Reference Remus, Mathis and Zahn2012) in their construction of a tidal model for small and moderate spin rates based on hydrodynamical equations (see Fig. 4 in Remus et al. Reference Remus, Mathis and Zahn2012).

For gravitational tides the time-lag is then given by

(75)$$b_{\rm g} ({\rm \sigma} ) = k_2 {\rm \sigma} {\rm \Delta} t_{\rm g} ({\rm \sigma} ) = \displaystyle{{k_2} \over {2Q_n}} \left[ {\displaystyle{{\rm \sigma} \over n} + \left( {2\,{\rm sign}({\rm \sigma} ) - \displaystyle{{\rm \sigma} \over n}} \right){\rm \xi} ({\rm \sigma} )} \right]\;, $$

where ξ(σ) is a function varying between 0 and 1, used to make a smooth passage between the two regimes:

(76)$${\rm \xi} ({\rm \sigma} ) = \displaystyle{1 \over 2} + \displaystyle{1 \over 2}\tanh {\rm \;} \left( {\displaystyle{{|{\rm \sigma} |} \over n} - 2} \right).$$

In Fig. 6, we plot the above expressions as a function of the tidal frequency. The transition of regime occurs for σ_c=2n, and Q _n corresponds to the dissipation Q – factor during the constant phase.

Fig. 6. (a) ξ(σ) versus σ with σ_c/n=2; (b) behaviour of b _g(σ) using the interpolated model smoothed by ξ(σ).

For thermal atmospheric tidesFootnote ¹³, the time-lag Δt _a(σ) is simply obtained from expression (55):

(77)$${\rm \Delta} t_{\rm a} ({\rm \sigma} ) = {\rm \Delta} t_{\rm g} ({\rm \sigma} )\displaystyle{{16H_0 K_{\rm g} k_2} \over {K_{\rm a} F_{\rm s}}} {\rm \omega} _{\rm s}. $$

Choice of parameters

Since much is unknown about exoplanets, we need to make some assumptions in our simulations. As we already stated, we consider here V-type planets. We assume that these planets with masses lower than 10 M _⊕ have a structure similar to the terrestrial planets of the Solar System, and a dense atmosphere capable of influencing their spin evolution. Thus, for a given exoplanet, we adopt for the unknown parameters the corresponding value from VenusFootnote ¹⁴. Therefore, in our simulations we take the potential Love number k ₂=0.28, Q _n=50, the mean density, $\bar {{ \rm \rho}} = 5.24 \times 10^3 \;{\rm kg}\;{\rm m} ^{- 3} $, and the stellar energy that reaches the planet surface F _s=100 W m⁻². For the mantle dynamic ellipticity we adopt δE _d=1.3×10⁻⁵ and for the kinematic viscosity ν=10⁻⁶ m² s⁻¹. For the core radius R _c=3.1×10⁶ m, the planet structure coefficient C/mR ²=0.336, and the ratio between the core moment of inertia and moment of inertia C _c/C=0.084 (Yoder Reference Yoder1995, Reference Yoder1997).

In our simulations we start with P _in=1 or 2 days, a rotation period that is faster than the predicted present rotation. Since some of these planets orbit very close to the star, and hence end up with fast equilibrium rotations (e.g., planet 55 Cnc e, whose final rotation is estimated to be 0.733 day), we also tested the possibility of planets evolved from a slower rotation to a faster one, with P _in=25 days. The results of these simulations are presented in the next two sections.

Application to already known exoplanets

In Fig. 7, we show the spin evolution for the planets HD 40307 b and g with time. We adopted P _in=1 day and initial obliquities ranging from 0° to 180°, with a step of 5°. Although the initial rotation rate is the same for planets b and g, the initial ratio ω/n is much higher for planet g, because the mean motion n is smaller. As a consequence, the characteristic evolution timescale is about 100 yr for planet b, while it is 1 Gyr for planet g. Those values are in agreement with the ones computed using equation (48) (Table 1). After that time all trajectories reach a final equilibrium for the spin, but depending on the initial value of the obliquity, the exact time can be different. Contrary to what we could expect, initial obliquities close to 180°, usually take less time to reach the equilibrium than lower initial values. This can be understood, since for high obliquities there are more harmonics with significant value that contribute to the spin evolution (equations (10, 11)).

Fig. 7. Spin evolution with time for the planets HD 40307 b (left) and HD 40307 g (right), with P _in=1 day and initial obliquities ranging from ε=0° to 180°. We plot the obliquity (top) and ω/n (bottom) evolution. Each line represents a different initial obliquity value. The lower lines in the ω/n plot correspond to the initial obliquities closer to 180°.

At first glance it may seem that there are two final equilibrium rotation states: one for final obliquities equal to 0°, and another for 180°. However, all initial obliquities evolve into the 2π/ω₁⁺ final rotation state, the only possibility for these two planets (Table 1). For the initial obliquities close 0°, the rotation rate is always positive and the obliquity is decreased to zero degrees. For initial obliquities close 180° the obliquity evolves into this value, while the rotation rate slows down until zero and then increases in the other way, until it stabilizes at a negative value −2π/|ω₁⁺|. This is the same as having a planet with zero obliquity rotating prograde, that is, the couple (−ω, π−ε) is mathematically equivalent to (ω, ε). Therefore, there are two different paths, but they lead to the same final equilibrium. This has also been described for Venus (Correia & Laskar Reference Correia and Laskar2001, Reference Correia and Laskar2003b).

It is also interesting to note that for initial obliquities around 90°, the rotation rate decreases to a value lower than that of equilibrium, and then it increases again until equilibrium is reached. This behaviour is in agreement with expression (47).

Comparing the rotation rate and the obliquity evolutions shown in Fig. 7 we also note that the obliquity is initially slightly constant and takes more time to reach the equilibrium than the rotation rate. The explanation for this behaviour is given in the ‘Obliquity evolution’ subsection. As we can see from expression (41), the obliquity variation is inversely proportional to the rotation rate. Thus, for a fast initial rotation, the obliquity does not change much. This is particularly visible for HD 40307 g, since the initial rotation ω≫n. However, as the rotation rate decreases we observe a strong variation in the obliquity. Finally, the obliquity slowly decreases into zero, since dx/dt∝(1−x ²)⇔dε/dt∝−sinε (equation (40)).

In Fig. 8, we show the obliquity evolution with the rotation rate ω/n for many systems listed in Table 1 (including HD 40307 b and g), always starting with an initial rotational period P=1 day. We observe again the dichotomy in the obliquity evolution: it can only end up in 0° or 180°. This result was demonstrated in Correia et al. (Reference Correia, Laskar and Néron de Surgy2003) for Venus’ parameters, but it seems to remain valid for all Earth-sized planets, even for those in eccentric orbits. This is in agreement with the assumption ε≃0° that we made when looking for final states (subsection ‘Equilibrium final states for the rotation rate’).

Fig. 8. Obliquity evolution with the rotation rate for several Earth-sized planets taken from Table 1 with an initial rotation period of P _in=1 day: (a) to (f) HD 40307 b to g; (g) 55 Cnc e; (h) GJ 1214 b; (i) HD 215497 b; (j) μ Arae c; (k) GJ 667C c; (l) HD 85512 b.

As a result of the obliquity's dichotomy, it appears that all planets have two possible final states, one prograde and another retrograde. For most of them only one final rotation state is possible, we are just observing the equivalent couples (ω, ε) and (−ω, π−ε). However, in one case for HD 40307 f (Fig. 8(e)), we can indeed observe two final rotation possibilities, the states ω₁⁺ and ω₁⁻ (Table 1). As these two states are so close to each other, it is nonetheless very difficult to distinguish them in the figure.

For all planets shown in Fig. 8, we also observe that the obliquity evolves into its final position at 0° or 180° always before ω/n=0. The reason is that for very slow rotation rates, core–mantle friction becomes dominating over tidal effects (equation (42)), quickly driving the obliquity to its final position. Therefore, the transition from positive into negative rotation rates only occurs at obliquities very close to 0° or 180°.

In the case of 55 Cnc e (Fig. 8(g)), we initially have ω/n≈0.76, that is, the initial rotation period of the planet is slower than the orbital period. As a consequence, contrarily to the remaining examples in Fig. 8, the rotation rate needs to increase its value to reach the equilibrium value near the synchronous rotation. We indeed observe this behaviour for initial obliquities close to 0°. However, for initial obliquities higher than 30°, the rotation rate still decreases to a lower ratio ω/n<0.76, before evolving into to the final equilibrium rotation. The explanation for this behaviour can be found in equation (47): for non-zero obliquities, the equilibrium rotation rate is below the synchronous rotation, so as long as the obliquity remains high we have 0<ω/n<1.

In Fig. 9, we show again the obliquity evolution with the rotation rate ω/n for the same systems in Fig. 8, but now starting with a slower rotation period P _in=25 days. Since most planets discovered so far are close-in planets with orbital periods smaller than 25 days, we now have an initial rotation rates ω<n for these planets. As a consequence, we observe an identical behaviour as the one described for 55 Cnc e with P _in=1 day (Fig. 8(g)).

Fig. 9. Obliquity evolution with the rotation rate for several Earth-sized planets taken from Table 1 with an initial rotation period of P _in=25 days: (a) to (f) HD 40307 b to g; (g) 55 Cnc e; (h) GJ 1214 b; (i) HD 215497 b; (j) μArae c; (k) GJ 667C c; (l) HD 85512 b.

Although the spin evolution with ω<n (Fig. 9) is completely different from the situation with ω>n (Fig. 8), the final rotation states are exactly the same. Different initial conditions lead to alternative evolutionary paths, but the final picture is the same. However, when ω<n we observe that there are much more trajectories increasing its obliquity to 180°. This behaviour is related to the onset of core–mantle friction. From equation (43) we have that (1) this effect is stronger for slow rotation rates; (2) when this effect becomes dominating, the obliquity evolves into 180° if ε>90° (x _f=x ₀ω₀/ω_f). Indeed, in Fig. 9 we clearly observe that when ω/n~0 there is a sudden inversion of the obliquity trend, whose final evolution is dictated by the 90° threshold.

We conclude that gravitational tides control the initial stages of the evolution, and that thermal atmospheric tides are only important in the definition of the final equilibrium rotation. Core–mantle friction only becomes important in the slow rotation regime (ω~n), but for very slow rotations it rules over both tidal effects.

Earth-sized planets in the habitable zone

Since most existing terrestrial planets are very close to their parent stars and/or orbit low-mass stars, usually there is only one final equilibrium rotation rate, given by state ω₁⁺ (Fig. 5, Table 1). In the examples from Figs. 8 and 9, only for HD 40703 f two final states are possible, but yet they are very close to each other and hardly distinguishable. As a consequence, in the simulations from the previous section, we were unable to observe evolutionary paths leading to very different configurations, such as retrograde rotation.

In the ‘Generalization to Earth-sized planets’ subsection, we investigated the most favourable orbital parameters that lead to multi-final states configurations. We saw that this correspond to planets not very close to the star (equation (74)), which may coincide with the HZ of their systems. Therefore, in order to observe more interesting evolutionary behaviours, we have also performed simulations for hypothetical Earth-sized planets around a Sun-like star, starting with an orbital period of 2 days. The final results from the simulations are in agreement with the equilibria states listed in Table 2.

In Fig. 10, we plotted some representative examples for different eccentricity and semi-major axis values. Although the initial rotation period is the same in all simulations, the initial ratio ω/n increases with the semi-major axis (since the mean motion decreases as the planet is moved away from the star). We still observe that the only possible final obliquities are ε=0° and 180°.

Fig. 10. Obliquity evolution with the rotation rate for several hypothetical Earth-sized planets taken from Table 2 with an initial rotation period of P _in=2 days. The eccentricities are e=0.0 (left), e=0.1 (middle), and e=0.2 (right). The products semi-amjor axis times mass of the parent star are (from top to bottom) aM _*=0.10, aM _*=0.42, aM _*=0.50 and aM _*=0.60 (units in [AU M _⊙]).

For different semi-major axis the main feature is the number of trajectories with increasing or decreasing obliquity. For close-in planets (a<0.1 AU) the picture is very similar to the real Earth-sized planets studied in the ‘Application to already known exoplanets’ subsection, but as the planet moves away from the star (or the stellar mass increase), more obliquities evolve into 180°. This is a direct consequence of thermal atmospheric tides, that push the obliquity towards 180°, and whose relative strength with respect to the gravitational tides increases for distant planets around Sun-like stars.

For a specific semi-major axis value, the global picture does not vary much for different eccentricities. However, the final evolution can be quite different, since the number of final states and their specific values can change completely (Table 2). For planets in circular orbits (e=0.0) only two final states are possible, symmetrical around synchronous rotation (ω=n±ω_s) (Figs. 4 and 5). However, the rotation periods only become substantially different from the synchronous rotation when planets are sufficiently distant from the star, since the effect from thermal atmospheric tides is more significant there. In particular, one of these states becomes retrograde for aM _*>0.59 AU M _⊙ (equation (74)). Indeed, when aM _*=0.6 AU M _⊙ (Fig. 10(j)), we can already observe a retrograde rotation state, 2π/ω₂⁻=−821 days, obtained whenever the initial obliquity is higher than 60°. The detailed final evolution near ε≃0° is shown in Fig. 11(a). Note, however, that while for initial 60°⩽ε<120° the obliquity is reduced to zero and the final state is reached through negative rotation, for ε⩾120° the obliquity evolves into 180° and the rotation rate is stabilized at positive rotation of 821 days. These final evolution scenarios are very similar to the case of Venus (a=0.72 AU), for which the presently observed rotation is also retrograde (Correia & Laskar Reference Correia and Laskar2001, Reference Correia and Laskar2003b).

Fig. 11. Detail of the final obliquity evolution with the rotation rate for (a) a=0.60 AU, e=0.0 and (b) a=0.42 AU, e=0.2. Two and three final equilibrium rotation states can be distinguished, respectively.

For non-circular orbits, the final evolution of the planets becomes more interesting. For instance, when e=0.1, the simulation with a=0.42 AU (Fig. 10(e)) shows three final possibilities for the spin (all prograde): for initial obliquities ε⩽40° trajectories go to the final state 2π/ω₂⁺=64 days; for 40°<ε<100° they go to the final state 2π/ω₁⁺=68 days, and for ε⩾100° they go to the final state 2π/ω₁⁻=176 days.

Similarly, for e=0.2 and the same semi-major axis a=0.42 AU (Fig. 10(f)), we still have the same three final states, but with slightly different periods. In Fig. 11(b), we show the detailed final evolution near ε=0° in this case, where the three distinct equilibria can be observed. For e=0.2 with slightly higher semi-major axis or stellar mass (aM _*=0.5 AU M _⊙) (Fig. 10(i)), we still observe three final states, but the final state ω₁⁻ is replaced by ω₂⁻, while for aM _*=0.6 AU M _⊙ (Fig. 10(l)) we are back to two final states (the ω₁⁺ final state disappears).

When we look at hypothetical Earth-sized planets with three possible final states, we observe that their sequence is always the same: lower values of initial obliquities correspond to the ω₂⁺ state, intermediate initial values correspond to the ω₁⁺ state, and higher initial obliquities to the ω₁⁻ or ω₂⁻ final states. The final equilibrium of the planet depends on the value of the rotation rate ω/n when the obliquity reaches 0°. Indeed, in the final evolution zooms shown in Fig. 11 we observe that the rotation of the planet evolves into to the equilibrium state that is closer to the rotation rate when it reaches ε=0°.