© E. F. S. Valente and A. C. M. Correia 2022

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 super-Earths, 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 Earth-mass planets is also expected to be large, as shown by formation studies (e.g. Schlecker et al. 2021). The fact that Earth-like 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 M-dwarfs. 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, Earth-sized planets are most frequent around early M-dwarfs (Burn et al. 2021). Moreover, these stars are perfect for detecting low-mass 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 M-dwarfs have already been discovered, some of which are promising candidates to harbour life (e.g. Anglada-Escudé 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 M-dwarfs 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 longest-living stars in the Universe, which provides enough time for planetary and biological development (Shields et al. 2016). However, it also has some disadvantages. M-dwarfs 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 (semi-major axis and eccentricity) of a planet are very important for habit-ability, 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 semi-major 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 slow-rotating 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 fast-rotating 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 high-obliquity planets are hotter than their low-obliquity 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 M-dwarfs 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 multi-planet systems, the eccentricity is excited by mutual perturbations and can keep a non-zero value (Laskar et al. 2012; Barnes 2017). In multi-planet 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 multi-body systems, for Earth-like close-in planets in the habitable zone of M-dwarfs, non-zero equilibrium states can be easily broken and evolve to a near-zero 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 two-body problem and does not require perturbations from other companions.

In Sect. 2, we present a general model to study the tidal evolution of close-in 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 Earth-mass planet in the habitable zone of an M-dwarf 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) $$\begin{array}{ccc}{E}_{\text{orb}}=\frac{GMm}{2a},& \text{and}& H=\beta n{a}^2\sqrt{1-e^2}k\end{array},$$(1)

where G is the gravitational constant, a is the semi-major axis, e is the eccentricity, n = [G(M + m)a⁻³]^1/2 is the mean motion, ß = Mm(M + m)⁻¹ 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 point-mass 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 $$L=C\omega ,$$(2)

where C = ξmR² 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), $$V\left(r\right)=k_2\frac{GM}{R}{\left(\frac{R}{r}\right)}^3{\left(\frac{R}{r_{\star }}\right)}^3{P}_2\left(\widehat{r}\cdot {\widehat{r}}_{\star }\right),$$(3)

where r is the distance measured from the planet centre of mass, r = ||r|| is the norm, $\widehat{r}=r/r$ is the unit vector, r_⋆ is the position of the star, and P₂(x) = (3x² − 1)/2 is a Legendre polynomial. k₂ is the second Love number for potential; it is a complex number, and depends on the frequency of the perturbation, σ. We can decompose k₂ in its real and imaginary parts as $$k_2\left(\sigma \right)=a\left(\sigma \right)-i\text{b}\left(\sigma \right).$$(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 $$g\left(r\right)=-{\nabla }_{r}V\left(r\right).$$(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 $$\dot{H}=T=M{r}_{\star }\times g\left(r_{\star }\right),$$(6)

and, owing to the conservation of the total angular momentum, $$\dot{L}=\dot{H}=-T.$$(7)

The evolution of the orbital energy (power) is given by $${\dot{E}}_{\text{orb}}=M{\dot{r}}_{\star }\cdot g\left(r_{\star }\right).$$(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 $$\begin{array}{ccc}p=\frac{k\times s}{\sin \theta },& q=\frac{k-\cos \theta s}{\sin \theta },& \cos \theta =k\cdot s,\end{array}$$(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 $$\langle T\rangle =T_P\sin \theta p+T_q\sin {\theta }_q+T_{s}S,$$(10)

and for the average power $$\langle {\dot{E}}_{\text{orb}}\rangle =n{T}_E,$$(11)

with¹ $$\begin{array}{ll}{T}_s=\frac{3\mathcal{K}}{32}\hfill & {\displaystyle \sum _{k=-\infty }^{+\infty }\{2b\left(-kn\right)\left(1-x^2\right)\left[{\left(X_k^{-3,2}\right)}^2-{\left(X_k^{-3,-2}\right)}^2\right]}\hfill \\ \hfill & +2b\left(\omega -kn\right){\left[1+x\right]}^2\left(2-x\right){\left(X_k^{-3,2}\right)}^2\hfill \\ \hfill & -4x^3{\left(X_k^{-3,0}\right)}^2-{\left(1-x\right)}^2\left(2+x\right){\left(X_k^{-3,-2}\right)}^2]\hfill \\ \hfill & +b\left(2\omega -kn\right)[-4x\left(1-x^2\right){\left(X_k^{-3,-2}\right)}^2\hfill \\ \hfill & +{\left(1+x\right)}^3{\left(X_k^{-3,2}\right)}^2-{\left(1-x\right)}^3{\left(X_k^{-3,-2}\right)}^2]\},\hfill \end{array}$$(12) $$\begin{array}{ll}{T}_s=\frac{3\mathcal{K}}{32}\hfill & {\displaystyle \sum _{k=-\infty }^{+\infty }\{2b\left(\omega -kn\right)\left(1-x\right)[4x^2{\left(X_k^{-3,0}\right)}^2}\hfill \\ \hfill & +{\left(1+x\right)}^2{\left(X_k^{-3,2}\right)}^2+{\left(1-x\right)}^2{\left(X_k^{-3,-2}\right)}^2]\hfill \\ \hfill & +b\left(2\omega -kn\right)[4{\left(1-x^2\right)}^2{\left(X_k^{-3,0}\right)}^2\hfill \\ \hfill & +{\left(1+x\right)}^4{\left(X_k^{-3,2}\right)}^2+{\left(1-x\right)}^4{\left(X_k^{-3,-2}\right)}^2]\},\hfill \end{array}$$(13) $$\begin{array}{ll}{T}_E=\frac{\mathcal{K}}{64}\hfill & {\displaystyle \sum _{k=-\infty }^{+\infty }}k\{b\left(-kn\right)[\left(1-3x^2\right){\left(X_k^{-3,0}\right)}^2\hfill \\ \hfill & +9{\left(1-x^2\right)}^2\left({\left(X_k^{-3,-2}\right)}^2+{\left(X_k^{-3,2}\right)}^2\right)]\hfill \\ \hfill & +12b\left(\omega -kn\right)\left(1-x^2\right)[4x^2{\left(X_k^{-3,0}\right)}^2\hfill \\ \hfill & +{\left(1-x\right)}^2{\left(X_k^{-3,2}\right)}^2+{\left(1+x\right)}^2{\left(X_k^{-3,2}\right)}^2]\hfill \\ \hfill & +3b\left(2\omega -kn\right)[4\left(1-x^2\right){\left(X_k^{-3,0}\right)}^2\hfill \\ \hfill & +\left(1-x\right){\left(X_k^{-3,-2}\right)}^2+{\left(1+x\right)}^4{\left(X_k^{-3,2}\right)}^2]\},\hfill \end{array}$$(14)

where K = GM²R⁵a⁻⁶, x = cos θ, and $X_k^{-3,m}\left(e\right)$ 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²⁰ ≈ 10⁻¹⁴, which corresponds to the computer precision that we have in our simulations.

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 $$\frac{\dot{a}}{a}=\frac{2a}{GMm}\langle {\dot{E}}_{\text{orb}}\rangle =\frac{2\sqrt{1-e^2}}{H}{T}_E,$$(15) $$\begin{array}{ll}\dot{e}\hfill & =\frac{\left(1-e^2\right)}{2e}\frac{\dot{a}}{a}-\frac{\sqrt{1-e^2}}{\beta n{a}^{2}e}\langle T\rangle \cdot k\hfill \\ \hfill & =\frac{1-e^2}{He}\left(T_E\sqrt{1-e^2}-T_q{\sin }^2\theta -T_s\cos \theta \right),\hfill \end{array}$$(16)

with H = ||H||. Similarly, the secular evolution of the rotation is obtained from the rotational angular momentum (Eq. (2)) as $$\frac{\dot{\omega }}{\omega }=-\frac{1}{C\omega }\langle T\rangle \cdot s=-\frac{T_s}{L},$$(17)

with L = ||L||, while the secular evolution of the obliquity is obtained from cos θ = k s = H L/(HL) as $$\begin{array}{ll}\dot{\theta }\hfill & =\left[\left(\frac{1}{L}+\frac{\cos \theta }{H}\right)\langle T\rangle \cdot k-\left(\frac{1}{H}+\frac{\cos \theta }{L}\right)\langle T\rangle \cdot s\right]\csc \theta \hfill \\ \hfill & =\left[\left(\frac{1}{L}+\frac{\cos \theta }{H}\right)T_q-\frac{T_s}{H}\right]\sin \theta .\hfill \end{array}$$(18)

In general, we have L ≪ H (Eq. (23)), and therefore the previous expression for the obliquity can be simplified as $$\dot{\theta }\approx \frac{T_q}{L}\sin \theta .$$(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₂(σ) 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 short-time and the long-time 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 long-term 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) $$\begin{array}{ccc}{k}_2\left(\sigma \right)=k_{\text{f}}\frac{1+\text{i}\sigma {\tau }_e}{1+\text{i}\sigma \tau },& \text{with}& \tau ={\tau }_e+{\tau }_{{}^{\upsilon }},\end{array}$$(20)

where τ, τ_υ, and τ_e = η/µ are the total, the viscous, and the Maxwell relaxation times, respectively. Here, k_f = k₂(0) is the fluid Love number, which corresponds to the maximal deformation that is attained for long-term perturbations. Similarly, we can also define an elastic Love number, k_e = k₂(∞) = k_f τ_e/τ, which corresponds to the deformation for short-time 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)), $$b\left(\sigma \right)=k_{\text{f}}\frac{\sigma {\tau }_{\upsilon }}{1+{\left(\sigma \tau \right)}^2}=\left(k_{\text{f}}-k_e\right)\frac{\sigma \tau }{1+{\left(\sigma \tau \right)}^2}.$$(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²¹ 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 post-glacial rebound due to the last glaciation about 10⁴ 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 $$\begin{array}{ccc}\frac{\dot{a}}{a}~\dot{e}~\frac{\mathcal{K}}{H}& \ll & \frac{\dot{\omega }}{\omega }~\dot{\theta }~\frac{\mathcal{K}}{L},\end{array}$$(22)

since (Eqs. (1) and (2)) $$\frac{L}{H}=\frac{C\omega }{\beta n{a}^2\sqrt{1-e^2}}\approx \frac{\omega }{n}{\left(\frac{R}{a}\right)}^2\ll 1.$$(23)

As a result, the spin evolves much faster than the orbit and is allowed to reach metastable or pseudo-equilibrium states as long as the eccentricity is not completely zero. Therefore, in this section we assume a constant value for the semi-major 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; Ferraz-Mello 2013; Correia et al. 2014). This is equivalent to setting x = 1 in expression (13), and thus expression (17) simplifies as $$\frac{\dot{\omega }}{\omega }=-\frac{3\mathcal{K}}{2L}{\displaystyle \sum _{k=-\infty }^{+\infty }b\left(2\omega -kn\right){\left(X_k^{-3,2}\left(e\right)\right)}^2.}$$(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 $\dot{\omega }=0$ 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 close-in planets with small eccentricities. Similar figures for higher eccentricity values can be found in Correia et al. (2014) and Ferraz-Mello (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) $$\frac{\dot{\omega }}{\omega }\approx \frac{3\mathcal{K}}{L}\left(k_{\text{f}}-k_{\text{e}}\right)n\tau \left(f_1\left(e\right)\frac{\omega }{n}-f_2\left(e\right)\right),$$(25)

with ƒ₁(e) ≈ 1 + 15e²/2, and ƒ₂(e) ≈ 1 + 27e²/2. The exact equilibrium, obtained when $\dot{\omega }=0$, is then given by $$\frac{\omega }{n}=\frac{f_2\left(e\right)}{f_1\left(e\right)}\approx 1+6e^2,$$(26)

which is also known as the pseudo-equilibrium 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: $$\frac{\dot{\omega }}{\omega }\approx -\frac{3\mathcal{K}}{4L}\frac{\left(k_{\text{f}}-k_{\text{e}}\right)}{n\tau }{\displaystyle \sum _{k=-\infty }^{+\infty }\frac{{\left(X_k^{-3,2}\left(e\right)\right)}^2}{\omega /n-k/2}.}$$(27)

For all k/2 values, $\dot{\omega }$ undergoes some oscillation, and its amplitude depends on $X_k^{-3,2}\left(e\right)$. The rotation is locked whenever the oscillation amplitude surpasses the line $\dot{\omega }=0$. These new states occur at almost semi-integer 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 $X_k^{-3,2}\left(e\right)$ 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 $\dot{\omega }=0$ 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 $\dot{\omega }=0$, 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 $$\frac{\dot{\omega }}{\omega }=-\frac{3\mathcal{K}}{4L}b\left(2\omega \right)f_1\left(e\right).$$(28)

As a consequence, we always have $\dot{\omega }<0$, 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 non-zero obliquities, the possible equilibria for the rotation rate remain at the semi-integer values ω/n = k/2, because they correspond to the zeros of the b(σ) functions. However, as for the eccentricity, non-zero 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 $X_k^{-3,2}\left(e\right)$ are truncated at |k| ≤ 20.

Fig. 2 Variation of dω/dt with ω/n around the 5/2 resonance for nτ =103 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 $X_k^{-3,2}\left(e\right)$ 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 $X_k^{-3,2}\left(e\right)$ 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 $$\begin{array}{c}\dot{\theta }\approx \frac{3\mathcal{K}}{4L}{\displaystyle \sum _{k=-\infty }^{+\infty }\{b\left(\omega -kn\right)\left[{\left(X_k^{-3,2}\left(e\right)\right)}^2-{\left(X_k^{-3,0}\left(e\right)\right)}^2\right]}\\ +b\left(2\omega -kn\right){\left(X_k^{-3,2}\left(e\right)\right)}^2\}\delta \theta .\end{array}$$(29)

In Fig. 3, we plot the previous expression as a function of ω/n for different values of the eccentricity and nτ. For $\dot{\theta }\text{/}\delta \theta \text{<0}$, we have stable equilibrium at θ = 0, while for $\dot{\theta }\text{/}\delta \theta \text{>0}$, this point is unstable. For initial fast rotations (ω/n ≫ 1), we can simplify expression (29) as (Correia & Valente 2022) $$\dot{\theta }\approx \frac{3\mathcal{K}}{4L}b\left(2\omega \right)f_1\left(e\right)\delta \theta .$$(30)

We conclude that we always have $\dot{\omega }<0$, 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 $\dot{\theta }\text{/}\delta \theta =0$ 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) $$\dot{\theta }\approx \frac{3\mathcal{K}}{L}\left(k_{\text{f}}-k_{\text{e}}\right)n\tau \left(f_1\left(e\right)\frac{\omega }{2n}-f_2\left(e\right)\right)\delta \theta ,$$(31)

which gives for the exact position of the transition point, $$\frac{\omega }{n}=2\frac{f_2\left(e\right)}{f_1\left(e\right)}\approx 2\left(1+6e^2\right).$$(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 $\dot{\theta }=0$. 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 high-obliquity 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 high-obliquity 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⁵, the obliquity can therefore reach nearly 80^° in the 5/2 resonance and 65^° in the 7/2 resonance. These interesting high-obliquity stable equibria were previously unknown, and in principle can be attained for some evolutionary paths of the spin.

Fig. 4 High-obliquity 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 Tq = Ts = 0 (Eqs. (12) and (13)). The Hansen coefficients $X_k^{-3,m}\left(e\right)$ 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 non-zero obliquity equilibria.

In Fig. 5, we show some trajectories for the spin in the plane (ω/n, θ) for different values of nτ (between 10⁻² and 10⁵) and constant eccentricity² (we adopt 0.05 ≤ e ≤ 0.16, because close-in Earth-like planets around M-dwarfs 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 $\dot{\theta }<0$ are plotted in red (the obliquity is damped), while segments with $\dot{\theta }>0$ 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 pseudo-rotation 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²), 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³, 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³, 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 high-obliquity state (Fig. 4). In Fig. 5, these stable high-obliquity equilibria are observed for (ω/n ≃ 5/2, θ ≃ 65°−80°) when e = 0.16 and nτ = 10³, e = 0.12 and nτ = 10⁴, or e = 0.09 and nτ = 10⁵, and for (ω/n ≃ 7/2, θ ≃ 65°) when e = 0.16 and nτ = 10⁵.

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², 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 close-in planets that depend on the initial obliquity, eccentricity, and relaxation time. For instance, for e = 0.09 and nτ = 10⁴, 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 semi-major 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 spin-state of a close-in Earth-like 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 long-term evolution of the system. Each panel can be seen as a snapshot at a given eccentricity, evolving from the right-hand 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 low-order 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 105. This plot was obtained with initial ω/n = 3.93, and by integrating Eqs. (17) and (19), where the Hansen coefficients $X_k^{-3,m}\left(e\right)$ are truncated at |k| ≤ 20. Segments with $\dot{\theta }<0$ are plotted in red, while segments with $\dot{\theta }>0$ are plotted in blue. Fixed points are marked with a black dot.

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.86-day orbit around a ~5 Gyr-old M-dwarf 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 M-dwarfs, 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 Earth-like 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 semi-major axis does not change much over the course of the simulations, even for nτ = 10². Indeed, owing to the conservation of the orbital angular momentum (Eq. (1)), the final semi-major axis is given by α_f ≈ α(1 − e²), which results in a variation ofonly 4% for e = 0.2. In addition, the semi-majoraxis only impacts the evolution timescale through the parameters K and H (Eq. (22)). Therefore, the choice of the initial semi-major 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⁵. 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² to 10⁵.

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², which corresponds to the highest tidal dissipation, because $\dot{\omega }\propto {\tau }^{-1}$ (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³. We observe that, although the tidal dissipation is about ten times weaker than for the case with nτ = 10² (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² 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⁴. The evolution of the spin is similar to the previous case with nτ = 10³, 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 high-obliquity 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 high-order spin–orbit resonances (Fig. 1), which includes the 7/2 resonance when nτ = 10⁴. 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⁴.

In Fig. 9, we finally adopt nτ = 10⁵, which corresponds to the more realistic value of nτ for Earth-like 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⁵, 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⁵ are not shown, because the eccentricity is already nearly constant for nτ = 10⁵. 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. High-obliquity 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 Earth-like planets in the habitable zone of M-dwarf stars, such as Proximab (Anglada-Escudé 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τ = 102 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τ = 103 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τ = 104 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

Earth-like rocky planets in the habitable zone of M-dwarfs, 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 habitabil-ity 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: $$P_{\text{syn}}=\frac{2\pi }{\omega -n}=\frac{P_{\text{orb}}}{\omega /n-1}.$$(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 M-dwarfs. 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 pseudo-synchronous 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 two-body problem, tides damp the eccentricity, which in turn destabilises the different equilibria for the spin. However, low-mass planets in the habitable zone of M-dwarfs are usually found in multi-planet systems (e.g. Dressing & Charbonneau 2015). As a result, the eccentricity can be excited by mutual perturbations and keep a non-zero value despite the tidal dissipation (Laskar et al. 2012; Barnes 2017). In multi-planet systems, the obliquity can also be locked in a Cassini state where it retains a non-zero value (Colombo 1966). Strong tidal effects usually destabilise high-obliquity 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 radial-velocity 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 close-in 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 ELT-METIS integral field spectrograph (Quanz et al. 2015; Brandl et al. 2021), will combine high-contrast imaging with high-resolution spectroscopy, and allow us to determine the spin state of these planets.

Fig. 9 Secular evolution over time of Ross 128 b with nτ = 105 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/FIS-AST/7002/2020), PHOBOS (POCI-01-0145-FEDER-029932), and ENGAGE SKA (POCI-01-0145-FEDER-022217), 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

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 $X_k^{-3,2}\left(e\right)$ are truncated at |k| ≤ 20.
In the text

	Fig. 2 Variation of dω/dt with ω/n around the 5/2 resonance for nτ =103 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 $X_k^{-3,2}\left(e\right)$ 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 $X_k^{-3,2}\left(e\right)$ are truncated at |k| ≤ 20.
In the text

	Fig. 4 High-obliquity 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 Tq = Ts = 0 (Eqs. (12) and (13)). The Hansen coefficients $X_k^{-3,m}\left(e\right)$ 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 105. This plot was obtained with initial ω/n = 3.93, and by integrating Eqs. (17) and (19), where the Hansen coefficients $X_k^{-3,m}\left(e\right)$ are truncated at |k| ≤ 20. Segments with $\dot{\theta }<0$ are plotted in red, while segments with $\dot{\theta }>0$ 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τ = 102 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τ = 103 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τ = 104 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τ = 105 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