Issue 
A&A
Volume 665, September 2022



Article Number  A130  
Number of page(s)  12  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202244010  
Published online  21 September 2022 
Tidal excitation of the obliquity of Earthlike planets in the habitable zone of Mdwarf stars
^{1}
CFisUC, Departamento de Física, Universidade de Coimbra,
3004516
Coimbra, Portugal
email: acor@uc.pt
^{2}
IMCCE, Observatoire de Paris, PSL Université,
77 Av. DenfertRochereau,
75014
Paris, France
Received:
12
May
2022
Accepted:
23
June
2022
Closein planets undergo strong tidal interactions with the parent star that modify their spins and orbits. In the twobody problem, the final stage for tidal evolution is the synchronisation of the rotation and orbital periods, and the alignment of the planet spin axis with the normal to the orbit (zero planet obliquity). The orbital eccentricity is also damped to zero, but over a much longer timescale, that may exceed the lifetime of the system. For nonzero eccentricities, the rotation rate can be trapped in spin–orbit resonances that delay the evolution towards the synchronous state. Here we show that capture in some spin–orbit resonances may also excite the obliquity to high values rather than damp it to zero. Depending on the system parameters, obliquities of 60º–80º can be maintained throughout the entire lifetime of the planet. This unexpected behaviour is particularly important for Earthlike planets in the habitable zone of Mdwarf stars, as it may help to sustain temperate environments and thus more favourable conditions for life.
Key words: planetstar interactions / planets and satellites: dynamical evolution and stability / planets and satellites: terrestrial planets / celestial mechanics / astrobiology
© E. F. S. Valente and A. C. M. Correia 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model. Subscribe to A&A to support open access publication.
1 Introduction
The discovery of a large number of exoplanets over recent decades shows that they are common around main sequence dwarf stars (e.g. Winn 2018). The most abundant planets are the superEarths, that is, planets with greater masses than that of the Earth and smaller than Neptune (e.g. Schlichting 2018). However, there exists an observational bias towards large masses, and so the population of Earthmass planets is also expected to be large, as shown by formation studies (e.g. Schlecker et al. 2021). The fact that Earthlike planets may exist around other stars offers a great opportunity to explore the possibility of having life as we know it outside the Solar System.
A configuration like the Earth around the Sun is ideal when searching for habitable worlds. However, researchers have shown interest in studying systems where the host stars are smaller than our Sun, the Mdwarfs. These stars are the most abundant in the Milky Way, representing ~70% of all stars (Bochanski et al. 2010), and the occurrence rate of planets is higher than around other main sequence stars (Tuomi et al. 2019). Formation studies also show that temperate, Earthsized planets are most frequent around early Mdwarfs (Burn et al. 2021). Moreover, these stars are perfect for detecting lowmass planets via radial velocity and transit techniques (Tuomi et al. 2014; Shields et al. 2016) because of their smaller mass and size. Indeed, several planets in the habitable zone of Mdwarfs have already been discovered, some of which are promising candidates to harbour life (e.g. AngladaEscudé et al. 2013, 2016; Dittmann et al. 2017; Gillon et al. 2017; Bonfils et al. 2018; Reiners et al. 2018; Zechmeister et al. 2019).
Searching for habitable worlds around Mdwarfs has some advantages. Because of their lower luminosities, the region where planets can sustain liquid water, usually called the habitable zone (Kasting et al. 1993), is closer to the star. In addition, these are the longestliving stars in the Universe, which provides enough time for planetary and biological development (Shields et al. 2016). However, it also has some disadvantages. Mdwarfs stars are very active stars, especially in their earlier lifetime, which drives atmospheric escape (e.g. France et al. 2020). The only possibility to avoid this is the presence of a strong magnetic field on the planets (Lammer et al. 2007), or a late formation of secondary atmospheres (Kite & Barnett 2020). Another concern is that, as the habitable zone is closer to the star, the planets undergo strong tidal effects, which modify their spins and orbits.
The spin (rotation and obliquity) and the orbit (semimajor axis and eccentricity) of a planet are very important for habitability, because they control the heat distribution on the surface of the planet, known as insolation (e.g. Milankovitch 1941; Ward 1974). The stellar flux is inversely proportional to the square of the semimajor axis, and so this parameter is used to define the limits of the habitable zone. High eccentricities also change the stellar flux on average, because of the different distances at the pericentre and appocentre (e.g. Dobrovolskis 2007). Therefore, planets with eccentric orbits in the habitable zone may not be able to sustain liquid water on the surface for the whole orbital period (Bolmont et al. 2016).
The rotation rate of a planet affects the day–night cycle, but also the atmosphere through the Coriolis effect. For instance, for slowrotating planets (like Venus) this effect is weak, which promotes the creation of thick stationary clouds in the substellar region, reflecting the solar rays (e.g. Kumar Kopparapu et al. 2019). On the other hand, for fastrotating planets (like the Earth) this effect is stronger and the clouds are in bands on the tropics, which in turn reduces the albedo (e.g. Yang et al. 2014).
The obliquity controls the insolation distribution for a given latitude, and is therefore the main driver of the seasons. High obliquities promote a more balanced distribution of the stellar flux on the planet surface, and thus extend the size of the habitable areas. Indeed, for obliquities of around 50º and 130°, the equatorial and polar regions receive identical amounts of insolation on average, while for obliquities close to 90° the polar regions receive as much insolation as the equatorial regions when the obliquity is 0° (e.g. Dobrovolskis 2021). Moreover, global climate model simulations show that highobliquity planets are hotter than their lowobliquity counterparts because of ice–albedo feedbacks for cold climates, which extends the outer edge of the habitable zone (Armstrong et al. 2014; Colose et al. 2019).
The final state of tidal evolution corresponds to circular orbits, zero obliquity, and synchronous rotation (Hut 1980; Adams & Bloch 2015). In this configuration, the rotation period is equal to the orbital period, implying that one side of the planet always faces the star and becomes extremely hot, while the other side becomes very cold. Though not impossible, this extreme environment is not suitable for the development of life as we know it (Wordsworth 2015; Turbet et al. 2016). However, the eccentricity of planets in the habitable zone of Mdwarfs is expected to evolve slowly, which allows the emergence of spin–orbit resonances that delay the evolution to the synchronous state (Makarov et al. 2012; Correia et al. 2014). Moreover, for multiplanet systems, the eccentricity is excited by mutual perturbations and can keep a nonzero value (Laskar et al. 2012; Barnes 2017). In multiplanet systems, the obliquity can also be locked in a ‘Cassini state’, where it retains a nonzero value (Colombo 1966). These equilibria occur more often for low obliquities, but they can also reach high values for systems with companions of large mass or high mutual inclination (Correia 2015). Resonant spin–orbit coupling can also trigger large chaotic variations of the obliquity (Laskar & Robutel 1993; Su & Lai 2022a).
Although high obliquities are possible in multibody systems, for Earthlike closein planets in the habitable zone of Mdwarfs, nonzero equilibrium states can be easily broken and evolve to a nearzero obliquity value (Levrard et al. 2007; Millholland & Laughlin 2019; Su & Lai 2022b). In this paper we describe a new mechanism that can stabilise the obliquity of these planets at high values, which is valid in the twobody problem and does not require perturbations from other companions.
In Sect. 2, we present a general model to study the tidal evolution of closein planets obtained in Correia & Valente (2022). In Sect. 3, we analyse the secular evolution of the rotation and obliquity in a very general framework, and show under which circumstances the obliquity is allowed to grow to high values. In Sect. 4, we apply our model to an Earthmass planet in the habitable zone of an Mdwarf star, Ross 128 b (Bonfils et al. 2018), and confirm the predictions from our model. Finally, in the last section we summarise and discuss our results.
2 Model
We consider a system in a Keplerian orbit, composed of a star and a planet, with masses M and m, respectively. The orbital energy and angular momentum are respectively given by (e.g. Murray & Dermott 1999) (1)
where G is the gravitational constant, a is the semimajor axis, e is the eccentricity, n = [G(M + m)a^{−3}]^{1/2} is the mean motion, ß = Mm(M + m)^{−1} is the reduced mass, and k is the unit vector along the direction of H, which is normal to the orbit. The star is a pointmass object, while the planet is an extended body with radius R, that can be deformed under the action of tides. The planet rotates with angular velocity ω = ωs, where s is the unit vector along the direction of the spin axis. We assume that s is also the axis of maximal inertia (gyroscopic approximation), and so the rotational angular momentum is simply given by (2)
where C = ξmR^{2} is the principal moment of inertia of the planet and ξ an internal structure constant.
2.1 Tidal torque and power
Estimations of the spin and orbital evolution of the planet are based on a very general formulation of the tidal potential initiated by Darwin (1879). Tidal effects arise from differential and inelastic deformations of the planet due to the gravitational perturbations from the star. The distortion of the planet gives rise to a tidal potential (e.g. Lambeck 1980), (3)
where r is the distance measured from the planet centre of mass, r = r is the norm, is the unit vector, r_{⋆} is the position of the star, and P_{2}(x) = (3x^{2} − 1)/2 is a Legendre polynomial. k_{2} is the second Love number for potential; it is a complex number, and depends on the frequency of the perturbation, σ. We can decompose k_{2} in its real and imaginary parts as (4)
This partition is very useful when we write the equations of motion (see Sect. 2.2), because the imaginary part characterises the viscous phase lag of the material and is therefore directly related to the amount of energy dissipated by tides. In general, the deformation lags behind the perturbation and therefore the imaginary part is always negative, hence the minus sign in expression (4).
The tidal potential (Eq. (3)) creates a differential gravitational field around the planet given by (5)
The star itself interacts with this field; with mass M and located at r = r _{⋆}, it exerts a torque on the planet that modifies its spin and orbit. In an inertial frame we have (6)
and, owing to the conservation of the total angular momentum, (7)
The evolution of the orbital energy (power) is given by (8)
2.2 Secular equations of motion
In general, tidal effects slowly modify the spin and the orbit of the planet, in a timescale much longer than the orbital and precession periods of the system. We can then average the torque (Eq. (6)) and the power (Eq. (8)) over the mean anomaly and the argument of the pericentre, and obtain the equations of motion for the secular evolution of the system.
We let (p, q, s) be a cartesian reference frame, such that (9)
where p is aligned with the line of nodes between the equator of the planet and the orbital plane, and θ is the angle between these two planes, also known as the obliquity. Following Correia & Valente (2022), for the average torque, we therefore have (10)
and for the average power (11)
where K = GM^{2}R^{5}a^{−6}, x = cos θ, and are the Hansen coefficients, which depend only on the eccentricity (see Table 1). Throughout this work, we consider only terms with k ≤ 20, because for the maximal adopted eccentricity (e = 0.2), we get an error smaller than e^{20} ≈ 10^{−14}, which corresponds to the computer precision that we have in our simulations.
Hansen coefficients and up to e^{6}.
2.3 Orbital and spin evolution
The secular evolution of the orbital elements is obtained from the orbital energy and angular momentum (Eq. (1)) as (15) (16)
with H = H. Similarly, the secular evolution of the rotation is obtained from the rotational angular momentum (Eq. (2)) as (17)
with L = L, while the secular evolution of the obliquity is obtained from cos θ = k s = H L/(HL) as (18)
In general, we have L ≪ H (Eq. (23)), and therefore the previous expression for the obliquity can be simplified as (19)
2.4 Maxwell rheology
The tidal deformation is characterised by the Love number (Eq. (4)), which depends on the internal structure of the planet and is therefore subject to large uncertainties. For this reason, in order to compute k_{2}(σ) we need to adopt some rheological model for the deformation. A large variety of models exist, but the most commonly used are the constant−Q (e.g. Munk & MacDonald 1960), the linear model (e.g. Mignard 1979), the Maxwell model (e.g. Correia et al. 2014), and the Andrade model (e.g. Efroimsky 2012). Some models appear to be better adapted to certain situations, but there is no model that is globally employed. However, viscoelastic rheologies are usually accepted as more realistic for rocky planets, because they are able to simultaneously reproduce the shorttime and the longtime responses (e.g. Remus et al. 2012). A review of the main viscoelastic models can be found in Renaud & Henning (2018).
We adopt here a Maxwell viscoelastic model, which is particularly well suited to reproducing the longterm deformation of the planets. A material is called a Maxwell solid when it responds to stresses like a massless, damped harmonic oscillator (e.g. Turcotte & Schubert 2002). It is characterised by a rigidity µ (or shear modulus), and by an effective viscosity η. Over short timescales, the material behaves like an elastic solid, but over long periods of time it flows like a fluid. The Love number for the Maxwell model is given by (e.g. Darwin 1908) (20)
where τ, τ_{υ}, and τ_{e} = η/µ are the total, the viscous, and the Maxwell relaxation times, respectively. Here, k_{f} = k_{2}(0) is the fluid Love number, which corresponds to the maximal deformation that is attained for longterm perturbations. Similarly, we can also define an elastic Love number, k_{e} = k_{2}(∞) = k_{f} τ_{e}/τ, which corresponds to the deformation for shorttime perturbations. Both k_{f} and k_{e} depend only on the internal structure of the planet, and can usually be measured for Solar System rocky planets from the rotational and tidal deformations, respectively (e.g. Correia & Rodríguez 2013). Therefore (Eqs. (4) and (20)), (21)
For the Earth, we have k_{e} = 0.295 and k_{f} = 0.933 (Yoder 1995). The total relaxation time can be obtained as τ = τ_{e} k_{f} /k_{e} ≈ 3η/µ. The rigidity of the Earth’s mantle is relatively well constrained, µ ≈ 80 GPa (Karato & Wu 1993). However, the average effective viscosity is much more uncertain, η ~ 10^{21} Pas (Karato & Wu 1993), which gives τ_{e} ~ 400 yr. Indeed, the viscosity may vary by many orders of magnitude, depending on the depth and temperature (e.g. Kirby & Kronenberg 1987). In the case of the Earth, the surface postglacial rebound due to the last glaciation about 10^{4} yr ago is still ongoing, suggesting that the Earth’s mantle relaxation time is τ_{e} ≈ 4400 yr (Turcotte & Schubert 2002). On the other hand, for slightly warmer planets we may have lower values, such as τ_{e} ~ 50 yr (e.g. Makarov et al. 2012).
3 Spin dynamics
In this section, we analyse the dynamics of the system provided by the secular equations (Sect. 2.3). The main goal is to describe the possible equilibria and evolutionary paths for different tidal regimes (different τ values).
Regardless of the tidal model that we adopt, the final state of tidal evolution corresponds to circular orbits, zero obliquity, and synchronous rotation (Hut 1980; Adams & Bloch 2015). However, the evolution timescales for the orbit and spin are quite different. Indeed, from expressions (15)–(19), we have (22)
As a result, the spin evolves much faster than the orbit and is allowed to reach metastable or pseudoequilibrium states as long as the eccentricity is not completely zero. Therefore, in this section we assume a constant value for the semimajor axis and eccentricity, and solely study the spin evolution.
We limit our analysis to positive rotation rates (ω > 0), but we can easily extend it to negative rotations with the symmetry (ω, θ) ⇔ (−ω, π − θ). These two points do not correspond to the same physical state, but are equivalent from a dynamical point of view (Correia & Laskar 2001).
We also note that the secular evolution of the spin (Eqs. (17) and (19)) depends only on four parameters: the rotation ratio, ω/n, the obliquity, θ, the eccentricity, e, and the relative relaxation time, nτ (Eq. (21)). The amplitude K/L can change the evolution rate, but does not modify the spin dynamics.
3.1 Equilibrium rotation
The final evolution of the spin has zero obliquity (θ = 0), which is also an equilibrium point (Eq. (19)). Therefore, for simplicity previous studies usually fix θ = 0, that is, they adopt a planar model (e.g. Makarov et al. 2012; FerrazMello 2013; Correia et al. 2014). This is equivalent to setting x = 1 in expression (13), and thus expression (17) simplifies as (24)
In Fig. 1, we plot the above expression as a function of ω/n for different values of the eccentricity and nτ. The interception with the horizontal dotted line gives all the equilibria for the rotation rate. Negative slopes correspond to stable equilibria, where the rotation can remain locked for a given value of eccentricity. We adopt a maximal eccentricity e = 0.16, because in this study we are only interested in closein planets with small eccentricities. Similar figures for higher eccentricity values can be found in Correia et al. (2014) and FerrazMello (2015).
For nτ ≲ 1, we observe that a single stable equilibrium point exists, close to synchronous rotation. Indeed, for small tidal frequencies (nτ ≪ 1), the Maxwell model behaves as the linear model, and expression (24) becomes (e.g. Correia 2009) (25)
with ƒ_{1}(e) ≈ 1 + 15e^{2}/2, and ƒ_{2}(e) ≈ 1 + 27e^{2}/2. The exact equilibrium, obtained when , is then given by (26)
which is also known as the pseudoequilibrium rotation.
For nτ ≫ 1, we observe that additional stable equilibria appear. They cluster around values ω/n ≈ k/2, because expression (24) has extrema at these values: (27)
For all k/2 values, undergoes some oscillation, and its amplitude depends on . The rotation is locked whenever the oscillation amplitude surpasses the line . These new states occur at almost semiinteger values (k ∈ ℤ), and they are usually called spin–orbit resonances (e.g. Goldreich & Peale 1966; Correia & Delisle 2019). The number of these equilibria increases with nτ, but also with the eccentricity value, because the amplitude of the Hansen coefficient increases with e (Table 1). We note that, except for the synchronous resonance (k = 2), for the remaining stable equilibria with k > 2, we have for a rotation rate that is slightly smaller than the exact resonant value, that is, for ω = kn/2 − δω, with δω/n ≪1. Indeed, as the oscillation in the torque introduced by the synchronous resonance has the largest amplitude (see Table 1), the average of the oscillations in the resonances with k > 2 is shifted to a negative value. In Fig. 2, we show an example for the 5/2 resonance with different eccentricities. For e = 0.11, the torque is always negative, and therefore capture cannot occur in this resonance. For e ≥ 0.12, the torque surpasses the line , and the rotation can stabilise at the value with negative slope, which occurs for ω/n = 5/2 − δω/n.
Immediately after their formation, the planets are supposed to rotate fast (e.g. Kokubo & Ida 2007), that is, ω/n ≫ 1. Therefore, we can write b(2ω − kn) ≈ b(2ω) and expression (24) simplifies as (28)
As a consequence, we always have , which can also be seen in Fig. 1 for ω/n = 4. We thus approach the equilibria while decreasing from fast rotations (right hand side in Fig. 1), and higher order spin–orbit resonances are encountered first.
In Fig. 1, we observe that for the same nτvalue, stabilisation can occur at different spin–orbit resonances, depending on the eccentricity. We are assuming here that the eccentricity is constant, but in reality it also slowly evolves (Eq. (22)). As the eccentricity decreases, the higher order spin–orbit equilibria become unstable, and the rotation is sequentially captured in lower order resonances, until it reaches its final possibility for synchronous rotation (see Sect. 4). Therefore, the asynchronous equilibria are only metastable; however, depending on the system parameters, they can persist for a time longer than the age of the system.
The general expression for the equilibria of the rotation rate (T_{s} = 0) depends on the obliquity value (Eq. (13)). For nonzero obliquities, the possible equilibria for the rotation rate remain at the semiinteger values ω/n = k/2, because they correspond to the zeros of the b(σ) functions. However, as for the eccentricity, nonzero obliquities can also modify the amplitude of the torque and therefore lead to a different distribution of the available states for a given initial condition.
Fig. 1 Variation of dω/dt with ω/n for different values of nτ and four different eccentricities: e = 0.05 (red), e = 0.09 (green), e = 0.12 (blue), and e = 0.16 (orange). We adopt θ = 0, and the Hansen coefficients are truncated at k ≤ 20. 
Fig. 2 Variation of dω/dt with ω/n around the 5/2 resonance for nτ =10^{3} and four different eccentricities: e = 0.11 (blue), e = 0.12 (green), e = 0.13 (orange), and e = 0.14 (red). We adopt θ = 0, and the Hansen coefficients are truncated at k ≤ 20. 
Fig. 3 Variation of dθ/dt with ω/n for different values of nτ and four different eccentricities: e = 0.05 (red), e = 0.09 (green), e = 0.12 (blue), and e = 0.16 (orange). We adopt θ ≈ 0, and the Hansen coefficients are truncated at k ≤ 20. 
3.2 Equilibrium obliquity
Although θ = 0 is always an equilibrium point (Eq. (19)), it is not necessarily stable. For T_{q} > 0 (Eq. (12)), this point is unstable and the obliquity is allowed to grow. Therefore, in this section we linearise expression (19) for θ = δθ ≪ 1 in order to determine when the planar approximation is no longer valid. We get (29)
In Fig. 3, we plot the previous expression as a function of ω/n for different values of the eccentricity and nτ. For , we have stable equilibrium at θ = 0, while for , this point is unstable. For initial fast rotations (ω/n ≫ 1), we can simplify expression (29) as (Correia & Valente 2022) (30)
We conclude that we always have , which can also be seen in Fig. 3 for ω/n = 4. Therefore, θ = 0 is an unstable equilibrium point, and the obliquity increases to higher values. However, as the rotation slows down owing to tides, we may have for one or more ω/n values, which represents a transition in the stability of the point θ = 0.
For nτ ≲ 1, we observe that a single transition point exists, close to ω/n = 2. In the linear model approximation (nτ ≪ 1), expression (29) becomes (e.g. Correia 2009) (31)
which gives for the exact position of the transition point, (32)
This value is two times larger than the equilibrium rotation (Eq. (26)). As a result, although θ = 0 is initially unstable and the obliquity allowed to grow, it becomes stable as the rotation approaches its equilibrium value.
For nτ ≫ 1, we observe that multiple transition points are present, introduced by an oscillation of the torque around the spin–orbit resonances ω/n = k/2, corresponding to the zeros of the b(σ) functions. As a result, the equilibrium at θ = 0 toggles between unstable and stable, and the transition points almost coincide with the spin–orbit resonances (Fig. 1). In general, the torque on the obliquity (Eq. (29)) is negative for ω/n < k/2 and positive for ω/n > k/2 (Fig. 3). Therefore, for spin–orbit resonances, the obliquity θ = 0 is usually stable, because for k > 2, the equilibrium rotation occurs for ω/n = k/2 − δω/n < k/2 (Fig. 2). Nevertheless, there are a few exceptions, namely whenever the amplitude of the oscillation is not large enough to descend below the horizontal line with . In those cases, the torque is always positive, and θ = 0 remains unstable. In Fig. 3, we observe that this is always the case for ω/n = 7/2, and it also occurs for ω/n = 5/2 for some eccentricities and nτvalues. We therefore expect that the obliquity grows when the rotation rate is captured in these two specific spin–orbit resonances.
When the obliquity increases, expression (29) and Fig. 3 are no longer valid. θ = π is also an equilibrium point, because sin π = 0 (Eq. (19)), but additional equilibria at intermediate obliquities may exist. As the expressions of T_{q} (Eq. (12)) and T_{s} (Eq. (13)) are very sensitive to the value of ω/n, the additional equilibria for a given eccentricity and nτvalue can only be found numerically by determining the pair (ω/n, θ) that simultaneously satisfy T_{q} = T_{s} = 0. In Fig. 4, we compute the highobliquity stable equilibria as a function of the eccentricity for different nτvalues. We find that for ω/n ≈ 5/2 and ω/n ≈ 7/2, there exist highobliquity states. These states do not directly depend on the value of nτ, but they become unstable below a critical eccentricity, which depends on nτ. We observe that as we increase the nτvalue, the critical eccentricity becomes smaller. For nτ = 10^{5}, the obliquity can therefore reach nearly 80^{°} in the 5/2 resonance and 65^{°} in the 7/2 resonance. These interesting highobliquity stable equibria were previously unknown, and in principle can be attained for some evolutionary paths of the spin.
Fig. 4 Highobliquity stable equilibria as a function of eccentricity. The dots mark the critical eccentricity for nτ = 10^{α}, where the αvalue is given next to the dot. These states are found numerically by determining the pair (ω/n, θ) that simultaneously satisfy T_{q} = T_{s} = 0 (Eqs. (12) and (13)). The Hansen coefficients are truncated at k ≤ 20. 
3.3 Global spin evolution
The rotation rate and the obliquity evolution cannot be dissociated as they progress at a similar pace (Eq. (22)). Therefore, the exact evolution of the spin can only be obtained by simultaneously integrating Eqs. (17) and (19). This allows us to confirm the analytic predictions from Sects. 3.1 and 3.2, in particular, the existence of nonzero obliquity equilibria.
In Fig. 5, we show some trajectories for the spin in the plane (ω/n, θ) for different values of nτ (between 10^{−2} and 10^{5}) and constant eccentricity^{2} (we adopt 0.05 ≤ e ≤ 0.16, because closein Earthlike planets around Mdwarfs are usually observedwith small eccentricities). We startthe integrations with ω/n = 3.93 and different obliquity values. As tidal effects initially decrease the rotation rate (Eq. (28)), the evolutionary paths must be followed for decreasing values of ω/n, until they reach an equilibrium position. Final metastable states (i.e. fixed points for a constant eccentricity) are marked with a black dot. To better understand the obliquity behaviour, segments with are plotted in red (the obliquity is damped), while segments with are plotted in blue (the obliquity grows).
For nτ ≲ 1, the system is in the linear regime. For initial small obliquities and ω/n > 2, the obliquity increases (Eq. (30)). However, as the rotation rate decreases, θ = 0 becomes stable (Eq. (32)). There are no other stable equilibrium points for the obliquity in this regime, and therefore all trajectories evolve to zero obliquity in the end. The final rotation rate is given by the pseudorotation equilibrium (Eq. (26)).
For nτ ≫ 1, multiple resonant equilibria exist at ω/n = k/2 (k ∈ ℤ), even for almost zero eccentricity. In general, the number of possibilities increases with the eccentricity and nτ (Fig. 1), but also with the obliquity (Eq. (12)). The resonant entrapments are well visible by vertical lines, because although the rotation is stalled, the obliquity continues to evolve. Most of these lines are in red, meaning that the obliquity decreases for those resonances. In some of them, the obliquity is able to reach θ = 0, where the rotation rate stabilises (black dot). In other cases (e.g. for ω/n = 3 with nτ = 10^{2}), as the obliquity decreases, the resonance becomes unstable, and the rotation decreases into a lower order resonance.
As nτ increases, more resonant equilibria become available (Fig. 1). For nτ ≳ 10^{3}, the obliquity is able to grow at some resonances (blue vertical lines). This behaviour is observed for ω/n = 5/2 and 7/2, because the torque on the obliquity can be always positive for those resonances (see Sect. 3.2 and Fig. 3). In some cases, the high obliquity ends by destabilising the resonant equilibrium^{3}, and the rotation rate decreases to a lower order spin–orbit resonance. More interestingly, however, we confirm that it is possible to stabilise the spin in a highobliquity state (Fig. 4). In Fig. 5, these stable highobliquity equilibria are observed for (ω/n ≃ 5/2, θ ≃ 65°−80°) when e = 0.16 and nτ = 10^{3}, e = 0.12 and nτ = 10^{4}, or e = 0.09 and nτ = 10^{5}, and for (ω/n ≃ 7/2, θ ≃ 65°) when e = 0.16 and nτ = 10^{5}.
As we start the simulations with ω/n ≈ 4, most trajectories are initially captured in the nearby 7/2 and 3/1 spin–orbit resonances. Nevertheless, lower order resonances can still be attained when the 3/1 resonance becomes unstable for low obliquity. Lower order resonances can also be directly attained for initial obliquities higher than 150^{°}, although in this case the initial obliquity has to be finely adjusted. For nτ = 10^{2}, we observe that initial obliquities higher than 165° evolve into 180°, while the rotation rate becomes negative. These simulations are continued at θ = 0 with the rotation rate increasing into the synchronous state, because from a dynamical point of view, the pair (−ω, θ = π) is equivalent to the pair (ω, θ = 0).
In Fig. 5, we observe that there is a large variety of evolutionary scenarios for the spin of closein planets that depend on the initial obliquity, eccentricity, and relaxation time. For instance, for e = 0.09 and nτ = 10^{4}, trajectories with initial obliquity lower than 150° are at first captured in the 3/1 spin–orbit resonance, where the obliquity decreases until it reaches a value close to zero. At this point, the 3/1 resonance becomes unstable, and the rotation decreases. Subsequently, capture in the 5/2 resonance occurs, but here the obliquity increases. When the obliquity reaches about 40^{°}, the 5/2 resonance becomes unstable, and the rotation rate finally evolves into the 2/1 spin–orbit resonance, where it stabilises with zero obliquity.
The different evolution scenarios depicted in Fig. 5 are completely general, as they do not depend on the semimajor axis, the radius, or the masses. These parameters appear in K/L and thus only modify the spin evolution timescale (Eq. (22)). In theory, for a given nτvalue, it is therefore possible to predict the present spinstate of a closein Earthlike planet solely from knowledge of its current eccentricity.
The spin evolution is also affected by the orbital evolution, because the eccentricity is not constant (Eq. (16)). The eccentricity generally decreases over a longer timescale, and so Fig. 5 can still be used to predict the longterm evolution of the system. Each panel can be seen as a snapshot at a given eccentricity, evolving from the righthand side to the left. As a result, the stable black dots only correspond to metastable states, because they disappear for decreasing eccentricity. For small eccentricities only loworder spin–orbit resonances persist, and for zero eccentricity only the synchronous rotation is possible for all nτvalues (see Fig. 5 in Correia et al. 2014).
Fig. 5 Secular evolution of the spin in the plane (ω/n, θ) for different eccentricities and nτvalues. From left to right, the eccentricity increases from 0.05 to 0.16. From top to bottom, the product nτ increases from 10^{−2} to 10^{5}. This plot was obtained with initial ω/n = 3.93, and by integrating Eqs. (17) and (19), where the Hansen coefficients are truncated at k ≤ 20. Segments with are plotted in red, while segments with are plotted in blue. Fixed points are marked with a black dot. 
Parameters of the Ross 128 system (Bonfils et al. 2018).
4 Application to Ross 128 b
To get the complete spin and orbital evolution of a system over time, we need to attribute specific values to all orbital and physical parameters. One good candidate for habitability studies is Ross 128 b, an exoplanet with a minimal mass of m = 1.35 M_{⊗} in a 9.86day orbit around a ~5 Gyrold Mdwarf star with M = 0.168 M_{⊙} (Bonfils et al. 2018). Although the planet is much closer to its host star (a ≈ 0.05 au) than the Earth is to the Sun, because of the much smaller luminosity of Mdwarfs, it only receives about 1.4 times as much flux as the Earth from the Sun. Ross 128 b is therefore likely to be situated at the inner edge of the habitable zone (Bonfils et al. 2018; Souto et al. 2018). In this section, we adopt the parameters listed in Table 2, and follow the tidal evolution of this planet by integrating the full set of secular equations of motion (15)–(18) for 10 Gyr.
The initial spin state of Earthlike planets is unknown. A small number of large impacts at the end of their formation can modify the rotation period and the axis orientation (Dones & Tremaine 1993; Kokubo & Ida 2007). However, as planets contract from a larger protoplanet, in general it is expected that they rotate fast. For simplicity, here we assume in all simulations that the initial rotation period of Ross 128 b is 24 h, which gives an initial rotation rate of ω/n = 9.86. This value is not critical, because the rotation rate rapidly decreases and evolves into an equilibrium value (Sect. 3.1). We also assume the same initial obliquity θ = 10^{°} in all simulations, because obliquities smaller than 150° are expected to behave in a similar way (Fig. 5). Moreover, by adopting a small initial obliquity value, it is possible to make the obliquity growth stand out.
The semimajor axis does not change much over the course of the simulations, even for nτ = 10^{2}. Indeed, owing to the conservation of the orbital angular momentum (Eq. (1)), the final semimajor axis is given by α_{f} ≈ α(1 − e^{2}), which results in a variation ofonly 4% for e = 0.2. In addition, the semimajoraxis only impacts the evolution timescale through the parameters K and H (Eq. (22)). Therefore, the choice of the initial semimajor axis is not critical, and so for simplicity we adopt the present value as the initial one.
The present eccentricity of Ross 128 b is not very well constrained from the observations, e = 0.12 ± 0.09 (Bonfils et al. 2018), that is, it is compatible with zero, but it can also be higher than 0.2. However, this parameter is critical for the spin evolution of the planet (Fig. 5). Therefore, we run our simulations adopting different values for the initial eccentricity from 0.01 to 0.20 in order to explore all the different behaviours.
Finally, although Ross 128 b is slightly more massive than the Earth, we assume a similar rheology for tides. For the density, structure constant, and Love numbers, we adopt those from the Earth (Table 2). As discussed in Sect. 2.4, the relaxation time is poorly constrained, but a typical adopted value is τe ~ 400 yr, which gives nτ ~ 10^{5}. However, this nτvalue leads to almost no orbital evolution over gigayear timescales, and therefore in order to get a more complete view of the possible spin evolutions, we run our simulations using different values, from nτ = 10^{2} to 10^{5}.
The results of the simulations for the tidal evolution of Ross 128 b are plotted in Figs. 6–9, one for each value of nτ. We show the eccentricity (top), the ω/n ratio (middle), and the obliquity (bottom) as a function of time for different initial eccentricities. For reasons of clarity, we only show six trajectories per plot, which illustrate all the important observed evolutions.
In Fig. 6, we adopt nτ = 10^{2}, which corresponds to the highest tidal dissipation, because (Eq. (27)). As a result, the rotation rate reaches the spin–orbit resonance regime in less than 1 Myr, and the orbit circularises in about 5 Gyr. We zoom into the early evolution (<1 Myr), for a clearer view of the initial spin evolution (left hand side of Fig. 6). We observe that for initial e = 0.2, the rotation rate is captured in the 5/2 resonance, for initial e = 0.1, it is captured in the 3/2 resonance, while for the remaining eccentricities between these two, all trajectories were captured in the 2/1 resonance. As expected, higher initial eccentricities foster capture in higher order spin–orbit resonances (Fig. 1). Before resonance entrapment, the obliquity increases for all simulations (Eq. (30)). However, as soon as capture in resonance occurs, the obliquity drops quickly to zero, except for the 5/2 resonance, where it continues to increase. After the obliquity surpasses a certain value (~35°), the 5/2 resonance becomes unstable, and the rotation decreases to the 2/1 resonance, where the obliquity is brought to zero. We limit the final evolution to 0.5 Gyr (right hand side of Fig. 6), because as the eccentricity decreases, all spin–orbit resonances become unstable one by one, until the rotation finally synchronises. Beyond that date, the eccentricity decays exponentially to zero, while the rotation remains synchronous and the obliquity at zero degrees.
In Fig. 7, we adopt nτ = 10^{3}. We observe that, although the tidal dissipation is about ten times weaker than for the case with nτ = 10^{2} (Eq. (27)), the rotation rate evolution does not differ much. It is initially trapped in some spin–orbit resonance, but as the eccentricity decreases, they all become unstable one by one. All trajectories finally end in the synchronous resonance before 5 Gyr (the present estimated age of the star). After that point, the eccentricity decreases very slowly, remaining above zero for at least 50 Gyr. One difference with respect to the case with nτ = 10^{2} is that the resonance entrapment persists for a much longer time, because the eccentricity is also damped much more slowly. Another difference is that for the trajectories captured in the 5/2 resonance, the obliquity is able to grow to a steady threshold around 65° (Fig. 4). As the eccentricity decreases, the equilibrium obliquity slightly increases, until the 5/2 resonance is destabilised. At this stage, the rotation rate decreases until it is captured in the lower order 2/1 resonance and the obliquity damps to zero. We observe that the obliquity growth in the 5/2 resonance is very consistent. For initial eccentricities higher than 0.16, the rotation is initially captured in the 3/1 resonance and the obliquity becomes very close to zero. However, as the 3/1 resonance becomes unstable, the rotation rate is trapped in the 5/2 resonance, where the obliquity rises again and settles in the 65° threshold. The high obliquity state is maintained as long as the eccentricity is above its critical value e ≈ 0.15 for nτ = 103 (Fig. 4). 4
In Fig. 8, we adopt nτ = 10^{4}. The evolution of the spin is similar to the previous case with nτ = 10^{3}, but the evolution timescale is about ten times longer. The eccentricity therefore varies very slowly and does not change much over the 10 Gyr of the simulation. As a consequence, capture in each spin–orbit resonance can last for several gigayears and only trajectories with initial e < 0.1 evolve into the synchronous rotation before 10 Gyr. For trajectories captured in the 5/2 resonance, the obliquity grows to nearly 75° and can remain there up to 3 Gyr. These highobliquity states are very robust, as they correspond to fixed points and can only be destabilised when the eccentricity descends below the critical value of e ≈ 0.11 (Fig. 4). Large nτvalues permit capture in highorder spin–orbit resonances (Fig. 1), which includes the 7/2 resonance when nτ = 10^{4}. In this resonance, the obliquity is also allowed to grow to high values. In our simulations, this scenario is observed for initial e = 0.20, but the obliquity cannot stabilise, because the critical eccentricity for this resonance is e ≈ 0.21 for nτ = 10^{4}.
In Fig. 9, we finally adopt nτ = 10^{5}, which corresponds to the more realistic value of nτ for Earthlike planets (see Sect. 2.4). In this case, tidal dissipation is so weak that the initial eccentricity almost does not change over the 10 Gyr of the simulation. As a result, the evolution of the spin can be depicted from the general evolution map shown in the last row of Fig. 5. Indeed, most trajectories remain throughout the simulation in the spin– orbit resonance where they were initially captured, except for a few cases (e = 0.08 and e = 0.14). For the trajectories that were captured in the 5/2 and 7/2 resonance, the obliquity is allowed to grow and stabilise at 79^{°} and 66^{°}, respectively (Fig. 4). For nτ = 10^{5}, we also observe one capture in the higher order 4/1 spin–orbit resonance for e = 0.2. This trajectory is not covered by the possibilities shown in Fig. 5, but in the 4/1 resonance the obliquity also evolves to zero.
Evolution trajectories for nτ > 10^{5} are not shown, because the eccentricity is already nearly constant for nτ = 10^{5}. We expect that the number of spin–orbit resonances increases for larger values of nτ (beyond the 4/1 resonance), which remain stable throughout the entire lifetime of the system. Highobliquity states are also possible and persist for long timescales, as they require lower eccentricities to become unstable (Fig. 4).
In this section, we present the results for Ross 128 b, but we also run simulations for other Earthlike planets in the habitable zone of Mdwarf stars, such as Proximab (AngladaEscudé et al. 2016) and Teegarden’s c (Zechmeister et al. 2019). We observed similar results, which confirms that the general analysis described in Sect. 3 is robust and independent from a fine tuning of the system parameters.
Fig. 6 Secular evolution over time of Ross 128 b with nτ = 10^{2} and some initial eccentricity values between 0.10 and 0.20. We show the eccentricity (top), the ω/n ratio (middle) and the obliquity (bottom). In all simulations, the initial rotation period is 24 h (ω/n = 9.86), and the initial obliquity is θ = 10^{º}. We note the break in the timescale. 
Fig. 7 Secular evolution over time of Ross 128 b with nτ = 10^{3} and some initial eccentricity values between 0.10 and 0.20. We show the eccentricity (top), the ω/n ratio (middle) and the obliquity (bottom). In all simulations, the initial rotation period is 24 h (ω/n = 9.86), and the initial obliquity is θ = 10^{º}. We note the break in the timescale. 
Fig. 8 Secular evolution over time of Ross 128 b with nτ = 10^{4} and some initial eccentricity values between 0.05 and 0.20. We show the eccentricity (top), the ω/n ratio (middle) and the obliquity (bottom). In all simulations, the initial rotation period is 24 h (ω/n = 9.86), and the initial obliquity is θ = 10^{°}. We note the break in the timescale. 
5 Discussion and conclusions
Earthlike rocky planets in the habitable zone of Mdwarfs, such as Ross 128 b, are expected to present a viscoelastic rheology with relaxation times such that nτ ≫ 1. In this regime, a multitude of spin–orbit resonances arise, where the spin can be temporarily trapped. Large nτvalues increase the number of possible resonances, but also their lifetime, because the eccentricity damping timescale is longer. In addition, for the 5/2 and 7/2 resonances, the obliquity can grow to high values of between 60^{°} and 80^{°}, which can be maintained for several gigayears. These results are extremely important regarding the habitability of these planets, because the rotation period and the obliquity determine the insolation distribution on the surface of the planet (e.g. Milankovitch 1941). The rotation rate controls the length of the day–night cycle through the synodic day period: (33)
For synchronous rotation (ω/n = 1), we have P_{syn} = ∞, and so one side of the planet always faces the star and becomes extremely hot, while the other side becomes very cold. Therefore, synchronous rotation is usually evoked as one of the major obstacles for the habitability of planets around Mdwarfs. However, in a 2/1 resonance, we have P_{syn} = P_{orb}, which corresponds to P_{syn} = 9.86 day in the case of Ross 128 b, while for the 7/2 resonance, we get P_{syn} = 3.94 day. The obliquity impacts the insolation distribution for a given latitude. It is responsible for the seasons, and therefore, even for synchronous planets, high obliquities can also mimic day–night cycles with the duration of one year (Dobrovolskis 2009), which corresponds to 9.86 days in the case of Ross 128 b. The combination of asynchronous rotation with extreme obliquities may thus help to sustain temperate environments and more favourable conditions for life.
Our results were obtained with a very general formalism, and are therefore not restricted to the study of rocky planets. They can be extended to other tidal regimes (different nτvalues) or to other rheologies by simply modifying the Love number (Eq. (20)). In our analysis (Sect. 3), we already considered relaxation times such that nτ ≲ 1, which are expected for gaseous planets (e.g. Correia et al. 2014). In this regime, the Maxwell model behaves as the linear model, and the results are well known from previous studies (e.g. Correia & Laskar 2010): the rotation rate evolves to the pseudosynchronous value (Eq. (26)), while the obliquity evolves to zero. These results are less interesting from a habitability point of view, but can still be used to put constraints on global atmospheric circulation models for planets in this regime, such as hot Jupiters or warm Neptunes.
In the twobody problem, tides damp the eccentricity, which in turn destabilises the different equilibria for the spin. However, lowmass planets in the habitable zone of Mdwarfs are usually found in multiplanet systems (e.g. Dressing & Charbonneau 2015). As a result, the eccentricity can be excited by mutual perturbations and keep a nonzero value despite the tidal dissipation (Laskar et al. 2012; Barnes 2017). In multiplanet systems, the obliquity can also be locked in a Cassini state where it retains a nonzero value (Colombo 1966). Strong tidal effects usually destabilise highobliquity Cassini states (e.g. Levrard et al. 2007), but in the peculiar 5/2 and 7/2 resonances the tidal torque can help to sustain those equilibria.
The tidal equilibria for the spin (ω/n, θ) depends only on two key parameters: the eccentricity, e, and the relative relaxation time, τ. The eccentricity is a parameter that we can already access with the current observational techniques, in particular for planets detected with the radialvelocity method. τ is therefore the only totally unknown parameter, and needs to be explored in order to put constraints on the possible spin state of a closein planet. Inversely, if we are able to measure the spin of a given planet, we can put constraints on τ. Although not possible at present, a new generation of instruments, such as the ELTMETIS integral field spectrograph (Quanz et al. 2015; Brandl et al. 2021), will combine highcontrast imaging with highresolution spectroscopy, and allow us to determine the spin state of these planets.
Fig. 9 Secular evolution over time of Ross 128 b with nτ = 10^{5} and some initial eccentricity values between 0.05 and 0.20. We show the eccentricity (top), the ω/n ratio (middle) and the obliquity (bottom). In all simulations, the initial rotation period is 24 h (ω/n = 9.86), and the initial obliquity is θ = 10^{°}. We note the break in the timescale. 
Acknowledgements
We are grateful to the referee Anthony Dobrovolskis for helpful comments. This work was supported by CFisUC (UIDB/04564/2020 and UIDP/04564/2020), GRAVITY (PTDC/FISAST/7002/2020), PHOBOS (POCI010145FEDER029932), and ENGAGE SKA (POCI010145FEDER022217), funded by COMPETE 2020 and FCT, Portugal. The numerical simulations were performed at the Laboratory for Advanced Computing at University of Coimbra (https://www.uc.pt/lca).
References
 Adams, F. C., & Bloch, A. M. 2015, MNRAS, 446, 3676 [Google Scholar]
 AngladaEscudé, G., Tuomi, M., Gerlach, E., et al. 2013, A&A, 556, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 AngladaEscudé, G., Amado, P. J., Barnes, J., et al. 2016, Nature, 536, 437 [Google Scholar]
 Armstrong, J., Barnes, R., DomagalGoldman, S., et al. 2014, Astrobiology, 14, 277 [CrossRef] [Google Scholar]
 Barnes, R. 2017, Celest. Mech. Dyn. l Astron., 129, 509 [NASA ADS] [CrossRef] [Google Scholar]
 Bochanski, J. J., Hawley, S. L., Covey, K. R., et al. 2010, AJ, 139, 2679 [NASA ADS] [CrossRef] [Google Scholar]
 Bolmont, E., Libert, A.S., Leconte, J., & Selsis, F. 2016, A&A, 591, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bonfils, X., AstudilloDefru, N., Díaz, R., et al. 2018, A&A, 613, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Boué, G., Correia, A. C. M., & Laskar, J. 2016, Celest. Mech. Dyn. Astron., 126, 31 [CrossRef] [Google Scholar]
 Brandl, B., Bettonvil, F., van Boekel, R., et al. 2021, The Messenger, 182, 22 [NASA ADS] [Google Scholar]
 Burn, R., Schlecker, M., Mordasini, C., et al. 2021, A&A, 656, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Colombo, G. 1966, AJ, 71, 891 [Google Scholar]
 Colose, C. M., Genio, A. D. D., & Way, M. J. 2019, ApJ, 884, 138 [NASA ADS] [CrossRef] [Google Scholar]
 Correia, A. C. M. 2006, Earth Planet. Sci. Lett. , 252, 398 [CrossRef] [Google Scholar]
 Correia, A. C. M. 2009, ApJ, 704, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Correia, A. C. M. 2015, A&A, 582, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Correia, A. C. M., & Delisle, J.B. 2019, A&A, 630, A102 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Correia, A. C. M., & Laskar, J. 2001, Nature, 411, 767 [Google Scholar]
 Correia, A. C. M., & Laskar, J. 2010, Icarus, 205, 338 [NASA ADS] [CrossRef] [Google Scholar]
 Correia, A. C. M., & Rodríguez, A. 2013, ApJ, 767, 128 [NASA ADS] [CrossRef] [Google Scholar]
 Correia, A. C. M., & Valente, E. F. S. 2022, Celest. Mech. Dyn. Astron., 134, 24 [NASA ADS] [CrossRef] [Google Scholar]
 Correia, A. C. M., Boué, G., Laskar, J., & Rodríguez, A. 2014, A&A, 571, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Darwin, G. H. 1879, Phil. Trans. R. Soc. London, 170, 1 [CrossRef] [Google Scholar]
 Darwin, G. H. 1908, Scientific Papers (Cambridge: Cambridge University Press) [Google Scholar]
 Dittmann, J. A., Irwin, J. M., Charbonneau, D., et al. 2017, Nature, 544, 333 [NASA ADS] [CrossRef] [Google Scholar]
 Dobrovolskis, A. R. 2007, Icarus, 192, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Dobrovolskis, A. R. 2009, Icarus, 204, 1 [CrossRef] [Google Scholar]
 Dobrovolskis, A. R. 2021, Icarus, 363, 114297 [NASA ADS] [CrossRef] [Google Scholar]
 Dones, L., & Tremaine, S. 1993, Icarus, 103, 67 [NASA ADS] [CrossRef] [Google Scholar]
 Dressing, C. D., & Charbonneau, D. 2015, ApJ, 807, 45 [Google Scholar]
 Efroimsky, M. 2012, Celest. Mech. Dyn. Astron., 112, 283 [NASA ADS] [CrossRef] [Google Scholar]
 FerrazMello, S. 2013, Celest. Mech. Dyn. Astron., 116, 109 [Google Scholar]
 FerrazMello, S. 2015, Celest. Mech. Dyn. Astron., 122, 359 [Google Scholar]
 France, K., Duvvuri, G., Egan, H., et al. 2020, AJ, 160, 237 [Google Scholar]
 Gillon, M., Triaud, A. H. M. J., Demory, B.O., et al. 2017, Nature, 542, 456 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & Peale, S. 1966, AJ, 71, 425 [NASA ADS] [CrossRef] [Google Scholar]
 Hut, P. 1980, A&A, 92, 167 [NASA ADS] [Google Scholar]
 Karato, S.I., & Wu, P. 1993, Science, 260, 771 [NASA ADS] [CrossRef] [Google Scholar]
 Kasting, J. F., Whitmire, D. P., & Reynolds, R. T. 1993, Icarus, 101, 108 [Google Scholar]
 Kirby, S. H., & Kronenberg, A. K. 1987, Rev. Geophys., 25, 1219 [NASA ADS] [CrossRef] [Google Scholar]
 Kite, E. S., & Barnett, M. N. 2020, Proc. Natl. Acad. Sci., 117, 18264 [NASA ADS] [CrossRef] [Google Scholar]
 Kokubo, E., & Ida, S. 2007, ApJ, 671, 2082 [NASA ADS] [CrossRef] [Google Scholar]
 Kumar Kopparapu, R., Wolf, E. T., & Meadows, V. S. 2019, Characterizing Exoplanet Habitability (Tucson: University of Arizona Press) [Google Scholar]
 Lambeck, K. 1980, The Earth’s Variable Rotation: Geophysical Causes and Consequences (Cambridge: Cambridge University Press) [CrossRef] [Google Scholar]
 Lammer, H., Lichtenegger, H., Kulikov, Y., et al. 2007, Astrobiology, 7, 185 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J., & Robutel, P. 1993, Nature, 361, 608 [Google Scholar]
 Laskar, J., Boué, G., & Correia, A. C. M. 2012, A&A, 538, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Levrard, B., Correia, A. C. M., Chabrier, G., et al. 2007, A&A, 462, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Makarov, V. V., Berghea, C., & Efroimsky, M. 2012, ApJ, 761, 83 [NASA ADS] [CrossRef] [Google Scholar]
 Mignard, F. 1979, Moon Planets, 20, 301 [Google Scholar]
 Milankovitch, M. 1941, Kanon der Erdbestrahlung und seine Andwendung auf das Eiszeitenproblem (Serbian: Royal Serbian Academy) [Google Scholar]
 Millholland, S., & Laughlin, G. 2019, Nat. Astron., 3, 424 [Google Scholar]
 Munk, W. H., & MacDonald, G. J. F. 1960, The Rotation of the Earth; A Geophysical Discussion (Cambridge: Cambridge University Press) [Google Scholar]
 Murray, C. D., & Dermott, S. F. 1999, Solar System Dynamics (Cambridge: Cambridge University Press) [Google Scholar]
 Quanz, S. P., Crossfield, I., Meyer, M. R., Schmalzl, E., & Held, J. 2015, Int. J. Astrobiol., 14, 279 [Google Scholar]
 Reiners, A., Ribas, I., Zechmeister, M., et al. 2018, A&A, 609, L5 [EDP Sciences] [Google Scholar]
 Remus, F., Mathis, S., & Zahn, J.P. 2012, A&A, 544, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Renaud, J. P., & Henning, W. G. 2018, ApJ, 857, 98 [Google Scholar]
 Schlecker, M., Pham, D., Burn, R., et al. 2021, A&A, 656, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schlichting, H. E. 2018, in Handbook of Exoplanets, eds. H. J. Deeg, & J. A. Belmonte (Berlin: Springer), 141 [Google Scholar]
 Shields, A. L., Ballard, S., & Johnson, J. A. 2016, Phys. Rep., 663, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Souto, D., Unterborn, C. T., Smith, V. V., et al. 2018, ApJ, 860, L15 [Google Scholar]
 Su, Y., & Lai, D. 2022a, MNRAS, 513, 3302 [NASA ADS] [CrossRef] [Google Scholar]
 Su, Y., & Lai, D. 2022b, MNRAS, 509, 3301 [Google Scholar]
 Tuomi, M., Jones, H. R. A., Barnes, J. R., AngladaEscudé, G., & Jenkins, J. S. 2014, MNRAS, 441, 1545 [Google Scholar]
 Tuomi, M., Jones, H. R. A., Butler, R. P., et al. 2019, AAS, submitted [arXiv: 1906.04644] [Google Scholar]
 Turbet, M., Leconte, J., Selsis, F., et al. 2016, A&A, 596, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Turcotte, D. L., & Schubert, G. 2002, Geodynamics (Cambridge: Cambridge University Press) [CrossRef] [Google Scholar]
 Ward, W. R. 1974, J. Geophys. Res., 79, 3375 [NASA ADS] [CrossRef] [Google Scholar]
 Winn, J. N. 2018, in Handbook of Exoplanets, eds. H. J. Deeg, & J. A. Belmonte (Berlin: Springer), 195 [Google Scholar]
 Wordsworth, R. 2015, ApJ, 806, 180 [NASA ADS] [CrossRef] [Google Scholar]
 Yang, J., Boué, G., Fabrycky, D. C., & Abbot, D. S. 2014, ApJ, 787, L2 [NASA ADS] [CrossRef] [Google Scholar]
 Yoder, C. F. 1995, in Global Earth Physics: A Handbook of Physical Constants (Washington D.C: American Geophysical Union), 1 [Google Scholar]
 Zechmeister, M., Dreizler, S., Ribas, I., et al. 2019, A&A, 627, A49 [NASA ADS] [EDP Sciences] [Google Scholar]
We do not provide the expression for T_{p}, because in our case this component changes neither the norm of the angular momentum vectors nor the angle between them (Sect. 2.3). This is also why we do not take into account the effect of the rotational deformation, as it only contributes to the torque with a component along p (e.g. Correia 2006).
Figure 5 in this paper is similar to Fig. 2 in Boué et al. (2016), who adopted nτ ≤ 10^{2} and different eccentricities (e = 0.0, 0.3 and 0.6). These choices did not allow us to spot some interesting behaviours, such as the stable nonzero obliquity equilibrium states.
A shortterm small obliquity excitation in the 7/2 resonance was originally reported by Boué et al. (2016) while using nτ = 10^{2} and e = 0.3. The obliquity climbs during about 15 Myr up to θ ≈ 15°, after which the 7/2 resonance becomes unstable.
All Tables
All Figures
Fig. 1 Variation of dω/dt with ω/n for different values of nτ and four different eccentricities: e = 0.05 (red), e = 0.09 (green), e = 0.12 (blue), and e = 0.16 (orange). We adopt θ = 0, and the Hansen coefficients are truncated at k ≤ 20. 

In the text 
Fig. 2 Variation of dω/dt with ω/n around the 5/2 resonance for nτ =10^{3} and four different eccentricities: e = 0.11 (blue), e = 0.12 (green), e = 0.13 (orange), and e = 0.14 (red). We adopt θ = 0, and the Hansen coefficients are truncated at k ≤ 20. 

In the text 
Fig. 3 Variation of dθ/dt with ω/n for different values of nτ and four different eccentricities: e = 0.05 (red), e = 0.09 (green), e = 0.12 (blue), and e = 0.16 (orange). We adopt θ ≈ 0, and the Hansen coefficients are truncated at k ≤ 20. 

In the text 
Fig. 4 Highobliquity stable equilibria as a function of eccentricity. The dots mark the critical eccentricity for nτ = 10^{α}, where the αvalue is given next to the dot. These states are found numerically by determining the pair (ω/n, θ) that simultaneously satisfy T_{q} = T_{s} = 0 (Eqs. (12) and (13)). The Hansen coefficients are truncated at k ≤ 20. 

In the text 
Fig. 5 Secular evolution of the spin in the plane (ω/n, θ) for different eccentricities and nτvalues. From left to right, the eccentricity increases from 0.05 to 0.16. From top to bottom, the product nτ increases from 10^{−2} to 10^{5}. This plot was obtained with initial ω/n = 3.93, and by integrating Eqs. (17) and (19), where the Hansen coefficients are truncated at k ≤ 20. Segments with are plotted in red, while segments with are plotted in blue. Fixed points are marked with a black dot. 

In the text 
Fig. 6 Secular evolution over time of Ross 128 b with nτ = 10^{2} and some initial eccentricity values between 0.10 and 0.20. We show the eccentricity (top), the ω/n ratio (middle) and the obliquity (bottom). In all simulations, the initial rotation period is 24 h (ω/n = 9.86), and the initial obliquity is θ = 10^{º}. We note the break in the timescale. 

In the text 
Fig. 7 Secular evolution over time of Ross 128 b with nτ = 10^{3} and some initial eccentricity values between 0.10 and 0.20. We show the eccentricity (top), the ω/n ratio (middle) and the obliquity (bottom). In all simulations, the initial rotation period is 24 h (ω/n = 9.86), and the initial obliquity is θ = 10^{º}. We note the break in the timescale. 

In the text 
Fig. 8 Secular evolution over time of Ross 128 b with nτ = 10^{4} and some initial eccentricity values between 0.05 and 0.20. We show the eccentricity (top), the ω/n ratio (middle) and the obliquity (bottom). In all simulations, the initial rotation period is 24 h (ω/n = 9.86), and the initial obliquity is θ = 10^{°}. We note the break in the timescale. 

In the text 
Fig. 9 Secular evolution over time of Ross 128 b with nτ = 10^{5} and some initial eccentricity values between 0.05 and 0.20. We show the eccentricity (top), the ω/n ratio (middle) and the obliquity (bottom). In all simulations, the initial rotation period is 24 h (ω/n = 9.86), and the initial obliquity is θ = 10^{°}. We note the break in the timescale. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.