# Tidal excitation of the obliquity of Earth-like planets in the habitable zone of M-dwarf stars

Ema F. S. Valente, Alexandre C. M. Correia

► To cite this version:

Ema F. S. Valente, Alexandre C. M. Correia. Tidal excitation of the obliquity of Earth-like planets in the habitable zone of M-dwarf stars. Astronomy and Astrophysics - A&A, 2022, 665, 10.1051/0004-6361/202244010. insu-03850268

HAL Id: insu-03850268

https://insu.hal.science/insu-03850268

Submitted on 13 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Tidal excitation of the obliquity of Earth-like planets in the habitable zone of M-dwarf stars

Ema F. S. Valente<fnref n="1" /> ID, and Alexandre C. M. Correia<fnref n="1" />,<fnref n="2" /> ID

Received 12 May 2022 / Accepted 23 June 2022

## ABSTRACT

Close-in planets undergo strong tidal interactions with the parent star that modify their spins and orbits. In the two-body 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 non-zero 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 Earth-like planets in the habitable zone of M-dwarf stars, as it may help to sustain temperate environments and thus more favourable conditions for life.

**Key words.** planet-star interactions – planets and satellites: dynamical evolution and stability – planets and satellites: terrestrial planets – celestial mechanics – astrobiology

## 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 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 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°,

<footnote n="1">¹ CFisUC, Departamento de Física, Universidade de Coimbra, 3004-516 Coimbra, Portugal
e-mail: acor@uc.pt</footnote>

<footnote n="2">² IMCCE, Observatoire de Paris, PSL Université, 77 Av. Denfert-Rochereau, 75014 Paris, France</footnote>
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 non-zero 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)

$$E_{\text{orb}} = -\frac{GMm}{2a}, \quad \text{and} \quad \mathbf{H} = \beta na^2 \sqrt{1 - e^2} \mathbf{k}, \qquad (1)$$

where $G$ is the gravitational constant, $a$ is the semi-major axis, $e$ is the eccentricity, $n = [G(M+m)a^{-3}]^{1/2}$ is the mean motion, $\beta = Mm(M+m)^{-1}$ is the reduced mass, and $\mathbf{k}$ is the unit vector along the direction of $\mathbf{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 $\omega = \omega 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, \qquad (2)$$

where $C = \xi m R^2$ is the principal moment of inertia of the planet and $\xi$ 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(\mathbf{r}) = -k_2 \frac{GM}{R} \left(\frac{R}{r}\right)^3 \left(\frac{R}{r_*}\right)^3 P_2(\hat{\mathbf{r}} \cdot \hat{\mathbf{r}}_*) \qquad (3)$$

where $\mathbf{r}$ is the distance measured from the planet centre of mass, $r = \|\mathbf{r}\|$ is the norm, $\hat{\mathbf{r}} = \mathbf{r}/r$ is the unit vector, $\mathbf{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, $\sigma$. We can decompose $k_2$ in its real and imaginary parts as

$$k_2(\sigma) = a(\sigma) - i b(\sigma). \qquad (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(\mathbf{r}) = -\nabla_r V(\mathbf{r}). \qquad (5)$$

The star itself interacts with this field; with mass $M$ and located at $\mathbf{r} = \mathbf{r}_*$, it exerts a torque on the planet that modifies its spin and orbit. In an inertial frame we have

$$\dot{\mathbf{H}} = \mathbf{T} = M \mathbf{r}_* \times g(\mathbf{r}_*), \qquad (6)$$

and, owing to the conservation of the total angular momentum,

$$\dot{\mathbf{L}} = -\dot{\mathbf{H}} = -\mathbf{T}. \qquad (7)$$

The evolution of the orbital energy (power) is given by

$$\dot{E}_{\text{orb}} = M \dot{\mathbf{r}}_* \cdot g(\mathbf{r}_*). \qquad (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

$$ \mathbf{p} = \frac{k \times s}{\sin \theta}, \quad \mathbf{q} = \frac{k - \cos \theta s}{\sin \theta}, \quad \cos \theta = \mathbf{k} \cdot \mathbf{s}, \qquad (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 \mathbf{T} \rangle = T_p \sin \theta \mathbf{p} + T_q \sin \theta \mathbf{q} + T_s \mathbf{s}, \qquad (10) $$

and for the average power

$$ \langle \dot{\mathbf{E}}_{\text{orb}} \rangle = n T_E, \qquad (11) $$

with<fnref n="1" />

$$ T_q = \frac{3K}{32} \sum_{k=-\infty}^{+\infty} \left\{ \begin{aligned} & 3b(-kn)(1-x^2)[(X_k^{-3,2})^2 - (X_k^{-3,-2})^2] \\ & + 2b(\omega - kn)[(1+x)^2(2-x)(X_k^{-3,2})^2 \\ & - 4x^3(X_k^{-3,0})^2 - (1-x)^2(2+x)(X_k^{-3,-2})^2] \\ & + b(2\omega - kn)[-4x(1-x^2)(X_k^{-3,0})^2 \\ & + (1+x)^3(X_k^{-3,2})^2 - (1-x)^3(X_k^{-3,-2})^2] \end{aligned} \right\}, \qquad (12) $$

$$ T_s = \frac{3K}{32} \sum_{k=-\infty}^{+\infty} \left\{ \begin{aligned} & 2b(\omega - kn)(1 - x^2)[4x^2(X_k^{-3,0})^2 \\ & + (1 + x)^2(X_k^{-3,2})^2 + (1 - x)^2(X_k^{-3,-2})^2] \\ & + b(2\omega - kn)[4(1 - x^2)^2(X_k^{-3,0})^2 \\ & + (1 + x)^4(X_k^{-3,2})^2 + (1 - x)^4(X_k^{-3,-2})^2] \end{aligned} \right\}, \qquad (13) $$

$$ T_E = \frac{K}{64} \sum_{k=-\infty}^{+\infty} k \left\{ \begin{aligned} & b(-kn) [4(1 - 3x^2)^2 (X_k^{-3,0})^2 \\ & + 9(1 - x^2)^2 [(X_k^{-3,-2})^2 + (X_k^{-3,2})^2]] \\ & + 12b(\omega - kn)(1 - x^2) [4x^2 (X_k^{-3,0})^2 \\ & + (1 - x)^2 (X_k^{-3,-2})^2 + (1 + x)^2 (X_k^{-3,2})^2] \\ & + 3b(2\omega - kn) [4(1 - x^2)^2 (X_k^{-3,0})^2 \\ & + (1 - x)^4 (X_k^{-3,-2})^2 + (1 + x)^4 (X_k^{-3,2})^2] \end{aligned} \right\}, \qquad (14) $$

<table><caption>**Table 1.** Hansen coefficients $X_k^{-3,0}(e)$ and $X_k^{-3,2}(e)$ up to $e^6$.</caption><thead><tr><th>k</th><th>X<sub>k</sub><sup>-3,0</sup>(e)</th><th>X<sub>k</sub><sup>-3,2</sup>(e)</th></tr></thead><tbody><tr><td>-6</td><td><sup>3167</sup>&frasl;<sub>320</sub>e<sup>6</sup></td><td>-</td></tr><tr><td>-5</td><td><sup>1773</sup>&frasl;<sub>256</sub>e<sup>5</sup></td><td>-</td></tr><tr><td>-4</td><td><sup>77</sup>&frasl;<sub>16</sub>e<sup>4</sup> + <sup>129</sup>&frasl;<sub>160</sub>e<sup>6</sup></td><td><sup>4</sup>&frasl;<sub>45</sub>e<sup>6</sup></td></tr><tr><td>-3</td><td><sup>53</sup>&frasl;<sub>16</sub>e<sup>3</sup> + <sup>393</sup>&frasl;<sub>256</sub>e<sup>5</sup></td><td><sup>81</sup>&frasl;<sub>1280</sub>e<sup>5</sup></td></tr><tr><td>-2</td><td><sup>9</sup>&frasl;<sub>4</sub>e<sup>2</sup> + <sup>7</sup>&frasl;<sub>4</sub>e<sup>4</sup> + <sup>141</sup>&frasl;<sub>64</sub>e<sup>6</sup></td><td><sup>1</sup>&frasl;<sub>24</sub>e<sup>4</sup> + <sup>7</sup>&frasl;<sub>240</sub>e<sup>6</sup></td></tr><tr><td>-1</td><td><sup>3</sup>&frasl;<sub>2</sub>e + <sup>27</sup>&frasl;<sub>16</sub>e<sup>3</sup> + <sup>261</sup>&frasl;<sub>128</sub>e<sup>5</sup></td><td><sup>1</sup>&frasl;<sub>48</sub>e<sup>3</sup> + <sup>11</sup>&frasl;<sub>768</sub>e<sup>5</sup></td></tr><tr><td>0</td><td>1 + <sup>3</sup>&frasl;<sub>2</sub>e<sup>2</sup> + <sup>15</sup>&frasl;<sub>8</sub>e<sup>4</sup> + <sup>35</sup>&frasl;<sub>16</sub>e<sup>6</sup></td><td>-</td></tr><tr><td>1</td><td><sup>3</sup>&frasl;<sub>2</sub>e + <sup>27</sup>&frasl;<sub>16</sub>e<sup>3</sup> + <sup>261</sup>&frasl;<sub>128</sub>e<sup>5</sup></td><td>-<sup>1</sup>&frasl;<sub>2</sub>e + <sup>1</sup>&frasl;<sub>16</sub>e<sup>3</sup> - <sup>5</sup>&frasl;<sub>384</sub>e<sup>5</sup></td></tr><tr><td>2</td><td><sup>9</sup>&frasl;<sub>4</sub>e<sup>2</sup> + <sup>7</sup>&frasl;<sub>4</sub>e<sup>4</sup> + <sup>141</sup>&frasl;<sub>64</sub>e<sup>6</sup></td><td>1 - <sup>5</sup>&frasl;<sub>2</sub>e<sup>2</sup> + <sup>13</sup>&frasl;<sub>16</sub>e<sup>4</sup> - <sup>35</sup>&frasl;<sub>288</sub>e<sup>6</sup></td></tr><tr><td>3</td><td><sup>53</sup>&frasl;<sub>16</sub>e<sup>3</sup> + <sup>393</sup>&frasl;<sub>256</sub>e<sup>5</sup></td><td><sup>7</sup>&frasl;<sub>2</sub>e - <sup>123</sup>&frasl;<sub>16</sub>e<sup>3</sup> + <sup>489</sup>&frasl;<sub>128</sub>e<sup>5</sup></td></tr><tr><td>4</td><td><sup>77</sup>&frasl;<sub>16</sub>e<sup>4</sup> + <sup>129</sup>&frasl;<sub>160</sub>e<sup>6</sup></td><td><sup>17</sup>&frasl;<sub>2</sub>e<sup>2</sup> - <sup>115</sup>&frasl;<sub>6</sub>e<sup>4</sup> + <sup>601</sup>&frasl;<sub>48</sub>e<sup>6</sup></td></tr><tr><td>5</td><td><sup>1773</sup>&frasl;<sub>256</sub>e<sup>5</sup></td><td><sup>845</sup>&frasl;<sub>48</sub>e<sup>3</sup> - <sup>32525</sup>&frasl;<sub>768</sub>e<sup>5</sup></td></tr><tr><td>6</td><td><sup>3167</sup>&frasl;<sub>320</sub>e<sup>6</sup></td><td><sup>533</sup>&frasl;<sub>16</sub>e<sup>4</sup> - <sup>13827</sup>&frasl;<sub>160</sub>e<sup>6</sup></td></tr><tr><td>7</td><td>-</td><td><sup>228347</sup>&frasl;<sub>3840</sub>e<sup>5</sup></td></tr><tr><td>8</td><td>-</td><td><sup>73369</sup>&frasl;<sub>720</sub>e<sup>6</sup></td></tr></tbody></table>

**Notes.** $X_k^{-3,-2}(e) = X_k^{-3,2}(e)$. The exact expression of these coefficients is given by $X_k^{-3,m}(e) = \pi^{-1} \int_0^\pi (a/r)^3 \exp(imv) \exp(-ikM) \, dM$.

where $K = GM^2R^5a^{-6}$, $x = \cos \theta$, and $X_k^{-3,m}(e)$ are the Hansen coefficients, which depend only on the eccentricity (see Table 1). Throughout this work, we consider only terms with $|k| \le 20$, because for the maximal adopted eccentricity ($e = 0.2$), we get an error smaller than $e^{20} \approx 10^{-14}$, 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{\mathbf{E}}_{\text{orb}} \rangle = \frac{2\sqrt{1-e^2}}{H} T_E, \qquad (15) $$

$$ \begin{aligned} \dot{e} &= \frac{(1 - e^2)}{2e} \frac{\dot{a}}{a} - \frac{\sqrt{1 - e^2}}{\beta na^2 e} \langle \mathbf{T} \rangle \cdot \mathbf{k} \\ &= \frac{1 - e^2}{He} (T_E \sqrt{1 - e^2} - T_q \sin^2 \theta - T_s \cos \theta), \end{aligned} \qquad (16) $$

with $H = \|\mathbf{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 \mathbf{T} \rangle \cdot \mathbf{s} = -\frac{T_s}{L}, \qquad (17) $$

with $L = \|\mathbf{L}\|$, while the secular evolution of the obliquity is obtained from $\cos \theta = \mathbf{k} \cdot \mathbf{s} = \mathbf{H} \cdot \mathbf{L}/(HL)$ as

$$ \begin{aligned} \dot{\theta} &= \left[ \left( \frac{1}{L} + \frac{\cos \theta}{H} \right) \langle \mathbf{T} \rangle \cdot \mathbf{k} - \left( \frac{1}{H} + \frac{\cos \theta}{L} \right) \langle \mathbf{T} \rangle \cdot \mathbf{s} \right] \csc \theta \\ &= \left[ \left( \frac{1}{L} + \frac{\cos \theta}{H} \right) T_q - \frac{T_s}{H} \right] \sin \theta. \end{aligned} \qquad (18) $$

<footnote n="1">¹ 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).</footnote>
In general, we have $L \ll H$ (Eq. (23)), and therefore the previous expression for the obliquity can be simplified as

$$ \dot{\theta} \approx \frac{T_q}{L} \sin \theta. \qquad (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(\sigma)$ 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 $\mu$ (or shear modulus), and by an effective viscosity $\eta$. 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)

$$ k_2(\sigma) = k_f \frac{1 + i\sigma\tau_e}{1 + i\sigma\tau}, \quad \text{with} \quad \tau = \tau_e + \tau_v, \qquad (20) $$

where $\tau$, $\tau_v$, and $\tau_e = \eta/\mu$ 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 long-term perturbations. Similarly, we can also define an elastic Love number, $k_e = k_2(\infty) = k_f \tau_e/\tau$, 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(\sigma) = k_f \frac{\sigma\tau_v}{1 + (\sigma\tau)^2} = (k_f - k_e) \frac{\sigma\tau}{1 + (\sigma\tau)^2}. \qquad (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 $\tau = \tau_e k_f/k_e \approx 3\eta/\mu$. The rigidity of the Earth's mantle is relatively well constrained, $\mu \approx 80$ GPa (Karato & Wu 1993). However, the average effective viscosity is much more uncertain, $\eta \sim 10^{21}$ Pa s (Karato & Wu 1993), which gives $\tau_e \sim 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^4$ yr ago is still ongoing, suggesting that the Earth's mantle relaxation time is $\tau_e \approx 4400$ yr (Turcotte & Schubert 2002). On the other hand, for slightly warmer planets we may have lower values, such as $\tau_e \sim 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 $\tau$ 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

$$ \frac{\dot{a}}{a} \sim \dot{e} \sim \frac{\mathcal{K}}{H} \ll \frac{\dot{\omega}}{\omega} \sim \dot{\theta} \sim \frac{\mathcal{K}}{L}, \qquad (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. \qquad (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 ($\omega > 0$), but we can easily extend it to negative rotations with the symmetry ($\omega, \theta \Leftrightarrow (-\omega, \pi - \theta)$). 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, $\omega/n$, the obliquity, $\theta$, the eccentricity, $e$, and the relative relaxation time, $n\tau$ (Eq. (21)). The amplitude $\mathcal{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 ($\theta = 0$), which is also an equilibrium point (Eq. (19)). Therefore, for simplicity previous studies usually fix $\theta = 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} \sum_{k=-\infty}^{+\infty} b(2\omega - kn) (X_k^{-3.2}(e))^2. \qquad (24) $$

In Fig. 1, we plot the above expression as a function of $\omega/n$ for different values of the eccentricity and $n\tau$. 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\tau \le 1$, we observe that a single stable equilibrium point exists, close to synchronous rotation. Indeed, for small tidal frequencies ($n\tau \ll 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} (k_f - k_e) n\tau \left( f_1(e) \frac{\omega}{n} - f_2(e) \right), \qquad (25) $$
<figure><img src="image_2.png" /><figcaption>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}(e) are truncated at |k| ≤ 20.</figcaption></figure>

with $f_1(e) \approx 1 + 15e^2/2$, and $f_2(e) \approx 1 + 27e^2/2$. The exact equilibrium, obtained when $\dot{\omega} = 0$, is then given by

$$ \frac{\omega}{n} = \frac{f_2(e)}{f_1(e)} \approx 1 + 6e^2, \qquad (26) $$

which is also known as the pseudo-equilibrium rotation.

For $n\tau \gg 1$, we observe that additional stable equilibria appear. They cluster around values $\omega/n \approx k/2$, because expression (24) has extrema at these values:

$$ \frac{\dot{\omega}}{\omega} \approx -\frac{3\mathcal{K}(k_f - k_c)}{4L n\tau} \sum_{k=-\infty}^{+\infty} \frac{(X_k^{-3,2}(e))^2}{\omega/n - k/2}. \qquad (27) $$

For all $k/2$ values, $\dot{\omega}$ undergoes some oscillation, and its amplitude depends on $X_k^{-3,2}(e)$. The rotation is locked whenever the oscillation amplitude surpasses the line $\dot{\omega} = 0$. These new states occur at almost semi-integer values ($k \in \mathbb{Z}$), and they are usually called spin-orbit resonances (e.g. Goldreich & Peale 1966; Correia & Delisle 2019). The number of these equilibria increases with $n\tau$, but also with the eccentricity value, because the amplitude of the Hansen coefficient $X_k^{-3,2}(e)$ 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 $\omega = kn/2 - |\delta\omega|$, with $|\delta\omega/n| \ll 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

<figure><img src="image_3.png" /><figcaption>Fig. 2. Variation of dω/dt with ω/n around the 5/2 resonance for nτ = 10³ 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}(e) are truncated at |k| ≤ 20.</figcaption></figure>

resonance. For $e \ge 0.12$, the torque surpasses the line $\dot{\omega} = 0$, and the rotation can stabilise at the value with negative slope, which occurs for $\omega/n = 5/2 - |\delta\omega/n|$.

Immediately after their formation, the planets are supposed to rotate fast (e.g. Kokubo & Ida 2007), that is, $\omega/n \gg 1$. Therefore, we can write $b(2\omega - kn) \approx b(2\omega)$ and expression (24) simplifies as

$$ \frac{\dot{\omega}}{\omega} = -\frac{3\mathcal{K}}{2L} b(2\omega)f_1(e). \qquad (28) $$

As a consequence, we always have $\dot{\omega} < 0$, which can also be seen in Fig. 1 for $\omega/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.
<figure><img src="image_4.png" /><figcaption>Fig. 3. Variation of $d\theta/dt$ with $\omega/n$ for different values of $n\tau$ and four different eccentricities: $e = 0.05$ (red), $e = 0.09$ (green), $e = 0.12$ (blue), and $e = 0.16$ (orange). We adopt $\theta \approx 0$, and the Hansen coefficients $X_k^{3,m}(e)$ are truncated at $|k| \le 20$.</figcaption></figure>

In Fig. 1, we observe that for the same $n\tau$-value, stabilization 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 $\omega/n = k/2$, because they correspond to the zeros of the $b(\sigma)$ 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.

## 3.2. Equilibrium obliquity

Although $\theta = 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 $\theta = \delta\theta \ll 1$ in order to determine when the planar approximation is no longer valid. We get

$$ \dot{\theta} \approx \frac{3\mathcal{K}}{4L} \sum_{k=-\infty}^{+\infty} \left\{ b(\omega - kn) \left[ (X_k^{-3,2}(e))^2 - (X_k^{-3,0}(e))^2 \right] + b(2\omega - kn) (X_k^{-3,2}(e))^2 \right\} \delta\theta. \quad (29) $$

In Fig. 3, we plot the previous expression as a function of $\omega/n$ for different values of the eccentricity and $n\tau$. For $\dot{\theta}/\delta\theta < 0$, we have stable equilibrium at $\theta = 0$, while for $\dot{\theta}/\delta\theta > 0$, this point is unstable. For initial fast rotations ($\omega/n \gg 1$), we can simplify expression (29) as (Correia & Valente 2022)

$$ \dot{\theta} \approx \frac{3\mathcal{K}}{4L} b(2\omega) f_1(e) \delta\theta. \quad (30) $$

We conclude that we always have $\dot{\theta}/\delta\theta > 0$, which can also be seen in Fig. 3 for $\omega/n = 4$. Therefore, $\theta = 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}/\delta\theta = 0$ for one or more $\omega/n$ values, which represents a transition in the stability of the point $\theta = 0$.

For $n\tau \ll 1$, we observe that a single transition point exists, close to $\omega/n = 2$. In the linear model approximation ($n\tau \ll 1$), expression (29) becomes (e.g. Correia 2009)

$$ \dot{\theta} \approx \frac{3\mathcal{K}}{L} (k_f - k_c) n\tau \left( f_1(e) \frac{\omega}{2n} - f_2(e) \right) \delta\theta, \quad (31) $$

which gives for the exact position of the transition point,

$$ \frac{\omega}{n} = 2 \frac{f_2(e)}{f_1(e)} \approx 2(1 + 6e^2). \quad (32) $$

This value is two times larger than the equilibrium rotation (Eq. (26)). As a result, although $\theta = 0$ is initially unstable and the obliquity allowed to grow, it becomes stable as the rotation approaches its equilibrium value.

For $n\tau \gg 1$, we observe that multiple transition points are present, introduced by an oscillation of the torque around the spin-orbit resonances $\omega/n = k/2$, corresponding to the zeros of
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** High-obliquity stable equilibria as a function of eccentricity. The dots mark the critical eccentricity for $n\tau = 10^4$, where the $\alpha$-value is given next to the dot. These states are found numerically by determining the pair $(\omega/n, \theta)$ that simultaneously satisfy $T_q = T_s = 0$ (Eqs. (12) and (13)). The Hansen coefficients $X_k^{-3m}(e)$ are truncated at $|k| \le 20$.</figcaption></figure>

the $b(\sigma)$ functions. As a result, the equilibrium at $\theta = 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 $\omega/n < k/2$ and positive for $\omega/n > k/2$ (Fig. 3). Therefore, for spin-orbit resonances, the obliquity $\theta = 0$ is usually stable, because for $k > 2$, the equilibrium rotation occurs for $\omega/n = k/2 - |\delta\omega/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 $\theta = 0$ remains unstable. In Fig. 3, we observe that this is always the case for $\omega/n = 7/2$, and it also occurs for $\omega/n = 5/2$ for some eccentricities and $n\tau$-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. $\theta = \pi$ is also an equilibrium point, because $\sin\pi = 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 $\omega/n$, the additional equilibria for a given eccentricity and $n\tau$-value can only be found numerically by determining the pair $(\omega/n, \theta)$ 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\tau$-values. We find that for $\omega/n \approx 5/2$ and $\omega/n \approx 7/2$, there exist high-obliquity states. These states do not directly depend on the value of $n\tau$, but they become unstable below a critical eccentricity, which depends on $n\tau$. We observe that as we increase the $n\tau$-value, the critical eccentricity becomes smaller. For $n\tau = 10^5$, the obliquity can therefore reach nearly 80° in the 5/2 resonance and 65° in the 7/2 resonance. These interesting high-obliquity stable equilibria were previously unknown, and in principle can be attained for some evolutionary paths of the spin.

### 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 $(\omega/n, \theta)$ for different values of $n\tau$ (between $10^{-2}$ and

$10^5$) and constant eccentricity<fnref n="2" /> (we adopt $0.05 \le e \le 0.16$, because close-in Earth-like planets around M-dwarfs are usually observed with small eccentricities). We start the integrations with $\omega/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 $\omega/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\tau \lesssim 1$, the system is in the linear regime. For initial small obliquities and $\omega/n > 2$, the obliquity increases (Eq. (30)). However, as the rotation rate decreases, $\theta = 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\tau \gg 1$, multiple resonant equilibria exist at $\omega/n = k/2$ ($k \in \mathbb{Z}$), even for almost zero eccentricity. In general, the number of possibilities increases with the eccentricity and $n\tau$ (Fig. 1), but also with the obliquity (Eq. (12)). The resonant entraps 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 $\theta = 0$, where the rotation rate stabilises (black dot). In other cases (e.g. for $\omega/n = 3$ with $n\tau = 10^2$), as the obliquity decreases, the resonance becomes unstable, and the rotation decreases into a lower order resonance.

As $n\tau$ increases, more resonant equilibria become available (Fig. 1). For $n\tau \gtrsim 10^3$, the obliquity is able to grow at some resonances (blue vertical lines). This behaviour is observed for $\omega/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<fnref n="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 high-obliquity state (Fig. 4). In Fig. 5, these stable high-obliquity equilibria are observed for $(\omega/n \approx 5/2, \theta \approx 65^\circ-80^\circ)$ when $e = 0.16$ and $n\tau = 10^3$, $e = 0.12$ and $n\tau = 10^4$, or $e = 0.09$ and $n\tau = 10^5$, and for $(\omega/n \approx 7/2, \theta \approx 65^\circ)$ when $e = 0.16$ and $n\tau = 10^5$.

As we start the simulations with $\omega/n \approx 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\tau = 10^2$, we observe that initial obliquities higher than 165° evolve into 180°, while the rotation rate becomes negative. These simulations are continued at $\theta = 0$ with the rotation rate increasing into the synchronous state, because from a dynamical point of view, the pair $(-\omega, \theta = \pi)$ is equivalent to the pair $(\omega, \theta = 0)$.

<footnote n="2">² Figure 5 in this paper is similar to Fig. 2 in Boué et al. (2016), who adopted $n\tau \le 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 non-zero obliquity equilibrium states.</footnote>

<footnote n="3">³ A short-term small obliquity excitation in the 7/2 resonance was originally reported by Boué et al. (2016) while using $n\tau = 10^2$ and $e = 0.3$. The obliquity climbs during about 15 Myr up to $\theta \approx 15^\circ$, after which the 7/2 resonance becomes unstable.</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Secular evolution of the spin in the plane $(\omega/n, \theta)$ for different eccentricities and $n\tau$-values. *From left to right*, the eccentricity increases from 0.05 to 0.16. *From top to bottom*, the product $n\tau$ increases from $10^{-2}$ to $10^5$. This plot was obtained with initial $\omega/n = 3.93$, and by integrating Eqs. (17) and (19), where the Hansen coefficients $X_k^{-3m}(e)$ are truncated at $|k| \le 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.</figcaption></figure>
<table><caption>**Table 2.** Parameters of the Ross 128 system (Bonfils et al. 2018).</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td><i>M</i></td><td>0.168 <i>M</i><sub>☉</sub></td></tr><tr><td><i>m</i></td><td>1.35 <i>M</i><sub>⊕</sub></td></tr><tr><td><i>a</i></td><td>0.0496 au</td></tr><tr><td>age</td><td>≥5 Gyr</td></tr><tr><td><i>R</i></td><td>1.12 <i>R</i><sub>⊕</sub></td></tr><tr><td><i>ξ</i></td><td>0.330</td></tr><tr><td><i>k</i><sub>e</sub></td><td>0.295</td></tr><tr><td><i>k</i><sub>f</sub></td><td>0.933</td></tr></tbody></table>

**Notes.** *R* is estimated from *m* assuming an Earth-like density. *ξ*, *k*<sub>e</sub> and *k*<sub>f</sub> are equal to the Earth's values (Yoder 1995).

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<sup>4</sup>, 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).

## 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*<sub>⊕</sub> in a 9.86-day orbit around a ~5 Gyr-old M-dwarf star with *M* = 0.168 *M*<sub>☉</sub> (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<sup>2</sup>. Indeed, owing to the conservation of the orbital angular momentum (Eq. (1)), the final semi-major axis is given by *a*<sub>f</sub> ≈ *a*(1 − *e*<sup>2</sup>), which results in a variation of only 4% for *e* = 0.2. In addition, the semi-major axis 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 τ<sub>e</sub> ~ 400 yr, which gives *n*τ ~ 10<sup>5</sup>. 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<sup>2</sup> to 10<sup>5</sup>.

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<sup>2</sup>, which corresponds to the highest tidal dissipation, because *ω* ~ τ<sup>-1</sup> (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
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Secular evolution over time of Ross 128 b with $n\tau = 10^2$ and some initial eccentricity values between 0.10 and 0.20. We show the eccentricity (top), the $\omega/n$ ratio (middle) and the obliquity (bottom). In all simulations, the initial rotation period is 24 h ($\omega/n = 9.86$), and the initial obliquity is $\theta = 10^\circ$. We note the break in the timescale.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Secular evolution over time of Ross 128 b with $n\tau = 10^3$ and some initial eccentricity values between 0.10 and 0.20. We show the eccentricity (top), the $\omega/n$ ratio (middle) and the obliquity (bottom). In all simulations, the initial rotation period is 24 h ($\omega/n = 9.86$), and the initial obliquity is $\theta = 10^\circ$. We note the break in the timescale.</figcaption></figure>

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\tau = 10^3$. We observe that, although the tidal dissipation is about ten times weaker than for the case with $n\tau = 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\tau = 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^\circ$ (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^\circ$ threshold. The high obliquity state is maintained as long as the eccentricity is above its critical value $e \approx 0.15$ for $n\tau = 10^3$ (Fig. 4).

In Fig. 8, we adopt $n\tau = 10^4$. The evolution of the spin is similar to the previous case with $n\tau = 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^\circ$ 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 \approx 0.11$ (Fig. 4). Large $n\tau$-values permit capture in high-order spin-orbit resonances (Fig. 1), which includes the 7/2 resonance when $n\tau = 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 \approx 0.21$ for $n\tau = 10^4$.

In Fig. 9, we finally adopt $n\tau = 10^5$, which corresponds to the more realistic value of $n\tau$ 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,
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Secular evolution over time of Ross 128 b with $n\tau = 10^4$ and some initial eccentricity values between 0.05 and 0.20. We show the eccentricity (top), the $\omega/n$ ratio (middle) and the obliquity (bottom). In all simulations, the initial rotation period is 24 h ($\omega/n = 9.86$), and the initial obliquity is $\theta = 10^\circ$. We note the break in the timescale.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Secular evolution over time of Ross 128 b with $n\tau = 10^5$ and some initial eccentricity values between 0.05 and 0.20. We show the eccentricity (top), the $\omega/n$ ratio (middle) and the obliquity (bottom). In all simulations, the initial rotation period is 24 h ($\omega/n = 9.86$), and the initial obliquity is $\theta = 10^\circ$. We note the break in the timescale.</figcaption></figure>

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\tau = 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\tau > 10^5$ are not shown, because the eccentricity is already nearly constant for $n\tau = 10^5$. We expect that the number of spin-orbit resonances increases for larger values of $n\tau$ (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 Proxima b (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.

## 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\tau \gg 1$. In this regime, a multitude of spin-orbit resonances arise, where the spin can be temporarily trapped. Large $n\tau$-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:

$$P_{\text{syn}} = \frac{2\pi}{\omega - n} = \frac{P_{\text{orb}}}{\omega/n - 1} \quad (33)$$

For synchronous rotation ($\omega/n = 1$), we have $P_{\text{syn}} = \infty$, 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_{\text{syn}} = P_{\text{orb}}$, which corresponds to $P_{\text{syn}} = 9.86$ day in the case of Ross 128 b, while for the 7/2 resonance, we get $P_{\text{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 rota-
tion 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 *nt*-values)
or to other rheologies by simply modifying the Love number
(Eq. (20)). In our analysis (Sect. 3), we already considered relax-
ation times such that *nt* ≲ 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 interest-
ing 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 dissipa-
tion (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 relax-
ation time, *τ*. The eccentricity is a parameter that we can already
access with the current observational techniques, in particu-
lar 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.

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

Adams, F. C., & Bloch, A. M. 2015, MNRAS, 446, 3676

Anglada-Escudé, G., Tuomi, M., Gerlach, E., et al. 2013, A&A, 556, A126

Anglada-Escudé, G., Amado, P. J., Barnes, J., et al. 2016, Nature, 536, 437

Armstrong, J., Barnes, R., Domagal-Goldman, S., et al. 2014, Astrobiology, 14, 277

Barnes, R. 2017, Celest. Mech. Dyn. l Astron., 129, 509

Bochanski, J. J., Hawley, S. L., Covey, K. R., et al. 2010, AJ, 139, 2679

Bolmont, E., Libert, A.-S., Leconte, J., & Selsis, F. 2016, A&A, 591, A106

Bonfils, X., Astudillo-Defru, N., Díaz, R., et al. 2018, A&A, 613, A25

Boué, G., Correia, A. C. M., & Laskar, J. 2016, Celest. Mech. Dyn. Astron., 126, 31

Brandl, B., Bettonvil, F., van Boekel, R., et al. 2021, The Messenger, 182, 22

Burn, R., Schlecker, M., Mordasini, C., et al. 2021, A&A, 656, A72

Colombo, G. 1966, AJ, 71, 891

Colose, C. M., Genio, A. D. D., & Way, M. J. 2019, ApJ, 884, 138

Correia, A. C. M. 2006, Earth Planet. Sci. Lett., 252, 398

Correia, A. C. M. 2009, ApJ, 704, L1

Correia, A. C. M. 2015, A&A, 582, A69

Correia, A. C. M., & Delisle, J.-B. 2019, A&A, 630, A102

Correia, A. C. M., & Laskar, J. 2001, Nature, 411, 767

Correia, A. C. M., & Laskar, J. 2010, Icarus, 205, 338

Correia, A. C. M., & Rodríguez, A. 2013, ApJ, 767, 128

Correia, A. C. M., & Valente, E. F. S. 2022, Celest. Mech. Dyn. Astron., 134, 24

Correia, A. C. M., Boué, G., Laskar, J., & Rodríguez, A. 2014, A&A, 571, A50

Darwin, G. H. 1879, Phil. Trans. R. Soc. London, 170, 1

Darwin, G. H. 1908, Scientific Papers (Cambridge: Cambridge University Press)

Dittmann, J. A., Irwin, J. M., Charbonneau, D., et al. 2017, Nature, 544, 333

Dobrovolskis, A. R. 2007, Icarus, 192, 1

Dobrovolskis, A. R. 2009, Icarus, 204, 1

Dobrovolskis, A. R. 2021, Icarus, 363, 114297

Dones, L., & Tremaine, S. 1993, Icarus, 103, 67

Dressing, C. D., & Charbonneau, D. 2015, ApJ, 807, 45

Efroimsky, M. 2012, Celest. Mech. Dyn. Astron., 112, 283

Ferraz-Mello, S. 2013, Celest. Mech. Dyn. Astron., 116, 109

Ferraz-Mello, S. 2015, Celest. Mech. Dyn. Astron., 122, 359

France, K., Duvvuri, G., Egan, H., et al. 2020, AJ, 160, 237

Gillon, M., Triaud, A. H. M. J., Demory, B.-O., et al. 2017, Nature, 542, 456

Goldreich, P., & Peale, S. 1966, AJ, 71, 425

Hut, P. 1980, A&A, 92, 167

Karato, S.-I., & Wu, P. 1993, Science, 260, 771

Kasting, J. F., Whitmire, D. P., & Reynolds, R. T. 1993, Icarus, 101, 108

Kirby, S. H., & Kronenberg, A. K. 1987, Rev. Geophys., 25, 1219

Kite, E. S., & Barnett, M. N. 2020, Proc. Natl. Acad. Sci., 117, 18264

Kokubo, E., & Ida, S. 2007, ApJ, 671, 2082

Kumar Kopparapu, R., Wolf, E. T., & Meadows, V. S. 2019, Characterizing Exoplanet Habitability (Tucson: University of Arizona Press)

Lambeck, K. 1980, The Earth's Variable Rotation: Geophysical Causes and Consequences (Cambridge: Cambridge University Press)

Lammer, H., Lichtenegger, H., Kulikov, Y., et al. 2007, Astrobiology, 7, 185

Laskar, J., & Robutel, P. 1993, Nature, 361, 608

Laskar, J., Boué, G., & Correia, A. C. M. 2012, A&A, 538, A105

Levrard, B., Correia, A. C. M., Chabrier, G., et al. 2007, A&A, 462, L5

Makarov, V. V., Berghea, C., & Efroimsky, M. 2012, ApJ, 761, 83

Mignard, F. 1979, Moon Planets, 20, 301

Milankovitch, M. 1941, Kanon der Erdbestrahlung und seine Andwendung auf das Eiszeitenproblem (Serbian: Royal Serbian Academy)

Millholland, S., & Laughlin, G. 2019, Nat. Astron., 3, 424

Munk, W. H., & MacDonald, G. J. F. 1960, The Rotation of the Earth; A Geophysical Discussion (Cambridge: Cambridge University Press)

Murray, C. D., & Dermott, S. F. 1999, Solar System Dynamics (Cambridge: Cambridge University Press)

Quanz, S. P., Crossfield, I., Meyer, M. R., Schmalzl, E., & Held, J. 2015, Int. J. Astrobiol., 14, 279

Reiners, A., Ribas, I., Zechmeister, M., et al. 2018, A&A, 609, L5

Remus, F., Mathis, S., & Zahn, J.-P. 2012, A&A, 544, A132

Renaud, J. P. & Henning, W. G. 2018, ApJ, 857, 98

Schlecker, M., Pham, D., Burn, R., et al. 2021, A&A, 656, A73

Schlichting, H. E. 2018, in Handbook of Exoplanets, eds. H. J. Deeg, & J. A. Belmonte (Berlin: Springer), 141

Shields, A. L., Ballard, S., & Johnson, J. A. 2016, Phys. Rep., 663, 1

Souto, D., Unterborn, C. T., Smith, V. V., et al. 2018, ApJ, 860, L15

Su, Y., & Lai, D. 2022a, MNRAS, 513, 3302

Su, Y., & Lai, D. 2022b, MNRAS, 509, 3301

Tuomi, M., Jones, H. R. A., Barnes, J. R., Anglada-Escudé, G., & Jenkins, J. S. 2014, MNRAS, 441, 1545

Tuomi, M., Jones, H. R. A., Butler, R. P., et al. 2019, AAS, submitted [arXiv:1906.04644]

Turbet, M., Leconte, J., Selsis, F., et al. 2016, A&A, 596, A112

Turcotte, D. L., & Schubert, G. 2002, Geodynamics (Cambridge: Cambridge University Press)

Ward, W. R. 1974, J. Geophys. Res., 79, 3375

Winn, J. N. 2018, in Handbook of Exoplanets, eds. H. J. Deeg, & J. A. Belmonte (Berlin: Springer), 195

Wordsworth, R. 2015, ApJ, 806, 180

Yang, J., Boué, G., Fabrycky, D. C., & Abbot, D. S. 2014, ApJ, 787, L2

Yoder, C. F. 1995, in Global Earth Physics: A Handbook of Physical Constants (Washington D.C: American Geophysical Union), 1

Zechmeister, M., Dreizler, S., Ribas, I., et al. 2019, A&A, 627, A49