# Tidal evolution of Earth-like planets in the habitable zone of low-mass stars

E. F. S. Valente, A. C. M. Correia, P. Auclair-Desrotour, M. Farhat, J. Laskar

► To cite this version:

E. F. S. Valente, A. C. M. Correia, P. Auclair-Desrotour, M. Farhat, J. Laskar. Tidal evolution of Earth-like planets in the habitable zone of low-mass stars. Astronomy & Astrophysics - A&A, 2024, 687, 10.1051/0004-6361/202348450. insu-04822204

HAL Id: insu-04822204

https://insu.hal.science/insu-04822204v1

Submitted on 6 Dec 2024

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 evolution of Earth-like planets in the habitable zone of low-mass stars

E. F. S. Valente<fnref n="1" /> Ⓓ, A. C. M. Correia<fnref n="1" />,<fnref n="2" /> Ⓓ, P. Auclair-Desrotour<fnref n="2" /> Ⓓ, M. Farhat<fnref n="2" /> Ⓓ, and J. Laskar<fnref n="2" /> Ⓓ

Received 31 October 2023 / Accepted 11 March 2024

## ABSTRACT

Earth-like planets in the habitable zone of low-mass stars undergo strong tidal effects that modify their spin states. These planets are expected to host dense atmospheres that can also play an important role in the spin evolution. On one hand, gravitational tides tend to synchronise the rotation with the orbital mean motion, but on the other hand, thermal atmospheric tides push the rotation away and may lead to asynchronous equilibria. Here, we investigate the complete tidal evolution of Earth-like planets by taking into account the effect of obliquity and eccentric orbits. We adopted an Andrade rheology for the gravitational tides and benchmarked the unknown parameters with the present rotation of Venus. We then applied our model to Earth-like planets, and we show that asynchronous rotation can be expected for planets orbiting stars with masses between 0.4 and 0.9 $M_⊕$ and semi-major axes between 0.2 and 0.7 au. Interestingly, we find that Earth-like planets in the habitable zone of stars with masses ~0.8 $M_⊕$ may end up with an equilibrium rotation of 24 h. We additionally find that these planets can also develop high obliquities, which may help sustain temperate environments.

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

## 1. Introduction

The number of detected low-mass exoplanets has grown quickly over the past decades. Nowadays, the most abundant planet classes are super-Earths and mini-Neptunes, that is, planets with masses larger than that of Earth and lower than that of Neptune (e.g. Fulton et al. 2017; Schlichting 2018). Because there is an observational bias towards large masses, a vast population of Earth-mass planets is also expected, as shown by formation studies (e.g. Schlecker et al. 2021). Many of these planets are currently found in the habitable zone (Kasting et al. 1993) of low-mass M-dwarf and K-dwarf stars (Burn et al. 2021). Therefore, we may wonder if they are similar to our Earth, and if so, which of these planets are possible candidates to harbour life as we know it.

When considering the habitability of a planet, not only does the host star have an influence but also the planetary properties, such as the composition of the atmosphere and the spin dynamics (e.g. Kopparapu et al. 2020). The mass of super-Earths is not large enough to allow these planets to retain a voluminous gaseous envelope, so they likely have a rocky composition. A thin atmospheric layer as it is observed on Earth and Venus is nevertheless expected (e.g. Demory et al. 2016; Komacek & Abbot 2019; Wordsworth & Kreidberg 2022), which can help sustain liquid water (Kasting et al. 1993; Kopparapu et al. 2013). The spin of these bodies controls the heat distribution on the surface from the incoming stellar radiation, known as insolation (e.g. Milankovitch 1941; Ward 1974), which strongly affects their surface temperature equilibrium and atmospheric general circulation (e.g. Forget & Leconte 2014; Dobrovolskis 2021). Therefore, the spin is a key quantity in predicting planet's climate.

Low-mass stars have a lower luminosity than Sun-like stars. As a result, the habitable zone is closer to the star, leading to significant star-planet interactions. Some of these stars can be very active, such as M-dwarfs in particular (Kopparapu et al. 2020), which can trigger significant atmospheric losses on close-in planets (e.g. Luger & Barnes 2015; Bolmont et al. 2017; France et al. 2020). Gravitational tides are also very strong, and they can modify the spin and the orbit of a planet until synchronous rotation is achieved (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<fnref n="1" />. Although synchronous rotation seems unfeasible for the development of life as we know it, it is not impossible (e.g. Yang et al. 2013; Turbet et al. 2016).

The habitable zone of K-dwarf stars is further away, so planets in such areas are more protected from strong stellar winds. However, the gravitational tidal interactions can still quickly modify their spins. Planets with dense atmospheres also experience thermal tides, which result from the heating of the atmosphere by the star. The torque induced by thermal atmospheric tides opposes that of gravitational tides, and thus it pushes the rotation away from the synchronisation. Although the mass of the atmosphere is often negligible with respect to the mass of the mantle, thermal tides can be of the same order of

<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, UMR8028 CNRS, Observatoire de Paris, PSL Université, 77 Av. Denfert-Rochereau, 75014 Paris, France</footnote>

<footnote n="1">¹ This is the case in the limit of thin atmospheres. However, the daytime temperature difference depends on the efficiency of the heat redistribution mechanisms, which has been examined in the literature using various modelling approaches (e.g. Wordsworth 2015; Koll & Abbot 2016; Auclair-Desrotour et al. 2022).</footnote>
magnitude (Dobrovolskis & Ingersoll 1980). Indeed, the peculiar retrograde asynchronous rotation of Venus is believed to result from the balance between gravitational and thermal atmospheric tides (Gold & Soter 1969; Ingersoll & Dobrovolskis 1978; Dobrovolskis 1980; Yoder 1997; Correia & Laskar 2001, 2003; Auclair-Desrotour et al. 2017b; Revol et al. 2023).

At present, thermal tides are negligible on Earth because they generate a torque of ~3.2 × 10¹⁵ J (e.g. Volland 1990), that is, less than 10% of the gravitational (ocean + solid) torque, a consequence of our planet being too distant from the Sun (Correia et al. 2008). Nevertheless, for Earth-like planets in the habitable zone of K-dwarf stars, thermal atmospheric tides can lead to asynchronous spin-orbit equilibria similar to that of Venus (Leconte et al. 2015). Moreover, thermal tides can also contribute to increasing the obliquity of these planets (e.g. Correia & Laskar 2003; Cunha et al. 2015). Moderate obliquities give rise to more temperate climates because the obliquity controls the insolation distribution for a given latitude and is therefore the main driver of the seasons (e.g. Dobrovolskis 2009, 2021).

In this paper, we revisit the spin evolution of Earth-like planets around low-mass stars mainly in the range of K-dwarf stars. For this purpose, we used a complete 3D dynamical model for the gravitational and thermal tides based on a vectorial formalism (Correia & Valente 2022; Valente & Correia 2023). We adopted the Andrade rheological model to handle the tidal deformation of the rocky mantle (e.g. Efroimsky 2012b; Gevorgyan et al. 2020) and the output of general circulation models (GCMs) to simulate the thermal response of the atmosphere (e.g. Leconte et al. 2015; Auclair-Desrotour et al. 2019); the models are believed to provide realistic descriptions of the atmosphere dynamics of rocky planets. Our aim is to better understand the possible end states for the spin of these planets and their impact on habitability as a function of stellar type, distance to the star, and orbital eccentricity.

In Sect. 2, we present the model that we used to study the tidal evolution of an Earth-mass planet around a single star. In Sect. 3, we apply this model to Venus in order to gauge the unknown parameters, assuming that Venus is currently on an equilibrium state (Correia & Laskar 2001). In Sect. 4, we apply our model to an Earth-twin planet (same mass and atmosphere properties of the Earth) and subject it to the environment of a low-mass star. We determine the equilibria for the spin for a wide range of stellar masses, semi-major axes, and eccentricities. In Sect. 5, we apply our model to the planet Kepler-1229 b, a super-Earth with radius 1.4 R⊕ in an 87-day orbit in the habitable zone of a star with mass 0.54 M⊕ (Morton et al. 2016), and compare with the general results from our model. Finally, in the last section, we summarise and discuss our results.

## 2. Model

We consider a system formed by a planet and a star moving in a Keplerian orbit, with masses $m$ and $M_★$, respectively. The orbital angular momentum and orbital energy are respectively given by (e.g. Murray & Dermott 1999)

$$G = \beta na^2 \sqrt{1 - e^2} k, \quad \text{and} \quad E = -\frac{GM_★m}{2a}, \qquad (1)$$

where $β = M_★m(M_★ + m)^{-1}$ is the reduced mass, $a$ is the semi-major axis, $e$ is the eccentricity, $n = [G(M_★ + m)a^{-3}]^{1/2}$ is the mean motion, $G$ is the gravitational constant and $k$ is the unit vector along the direction of $G$, which is normal to the orbit.

The star is a point-mass object with luminosity $L_★$, while the planet is an extended body composed by a mantle with radius $R$, and surrounded by a thin atmosphere. It rotates with angular velocity $ω = ωs$, where $s$ is the unit vector along the direction of the spin axis. We assumed that $s$ is also the axis of maximal inertia (gyroscopic approximation), so the rotational angular momentum is simply given by

$$L = Cω, \qquad (2)$$

where $C = ξmR^2$ is the principal moment of inertia of the planet and $ξ$ is an internal structure constant. We further assumed that the mantle of the planet and its atmosphere can be deformed under the action of gravitational and thermal tides, respectively.

### 2.1. Gravitational tides

Gravitational tides arise from differential and inelastic deformations of the mantle due to the gravitational perturbations of the star. The atmosphere also undergoes gravitational tides, but they can be neglected. The distortion of the planet gives rise to a tidal potential (e.g. Lambeck 1980),

$$V_g(\mathbf{r}) = -\hat{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 $r$ is the distance measured from the planet centre-of-mass, $r = ||r||$ is the norm, $\hat{\mathbf{r}} = \mathbf{r}/r$ is the unit vector, $\mathbf{r}_★$ is the position of the star, $P_2(x) = (3x^2 - 1)/2$ is a Legendre polynomial and $\hat{k}_2$ is the second Love number for potential.

The deformation of the mantle is characterised by the Love number, $\hat{k}_2$, which depends on the frequency of the perturbation, $\sigma$, and on the internal structure of the planet. In general, we can write (e.g. Lambeck 1980)

$$\hat{k}_2(\sigma) = \frac{k_f}{1 + \hat{\mu}(\sigma)}, \qquad (4)$$

where $k_f$ is the fluid Love number and $\hat{\mu}$ is an effective rigidity. For a homogenous body, $k_f = 3/2$, but more generally it depends on the mass distribution inside the planet and can be related to the structure constant $ξ$ through the Darwin-Radau equation (e.g. Jeffreys 1976). For Earth and Venus, we have $k_f = 0.933$ and $k_f = 0.928$, respectively (Yoder 1995).

In order to compute the effective rigidity, $\hat{\mu}(\sigma)$, we need to adopt some rheological model for the deformation. The exact response of rocky planets to stress is unknown due to the lack of information about their interiors (e.g. Mocquet et al. 2011). However, viscoelastic rheologies are usually accepted as more realistic, because they are able to simultaneously reproduce the short-time and the long-time responses (e.g. Remus et al. 2012; Efroimsky 2012a; Correia et al. 2014; Renaud & Henning 2018).

Transient creep of many materials at high temperatures obeys the Andrade's law (e.g. Louchet & Duval 2009). Moreover, geophysical data from body and surface waves, free oscillations, and Chandler wobble, which cover a large spectrum of frequencies, also suggest that the mantle dissipates energy according to an Andrade rheology (e.g. Andrade 1910; Anderson & Minster 1979; Efroimsky 2012a; Gevorgyan et al. 2020). Therefore, we adopt here the Andrade viscoelastic model. The effective rigidity is then given by (e.g. Efroimsky 2012a)

$$\hat{\mu}(\sigma) = \frac{\tau_v}{\tau_e} \left[ 1 - i(\sigma\tau_e)^{-1} + (i\sigma\tau_A)^{-\alpha}\Gamma(1 + \alpha) \right]^{-1}, \qquad (5)$$
where $\tau_v$ is the viscous relaxation time, $\tau_e = \eta/\mu$ is the Maxwell relaxation time, $\eta$ is the effective viscosity and $\mu$ is the rigidity. The parameter $\alpha$ is the time exponent of the Andrade transient creep strain law, which is an empirical adjustable constant that usually ranges from $0.2 \le \alpha \le 0.4$ for rocky planets (Castillo-Rogez et al. 2011). The quantity $\tau_A$ is the timescale associated with the Andrade creep and may be termed as the “Andrade” or the “anelastic” time. In general, we have $\tau_A \ge \tau_e$ (Makarov 2012), so for simplicity<fnref n="2" />, we assumed that $\tau_A = \tau_e$.

The effective rigidity $\hat{\mu}$ is a complex number (Eq. (5)), and hence also the Love number $\hat{k}_2$ (Eq. (4)). We can thus decompose $\hat{k}_2$ in its real and imaginary parts as

$$ \hat{k}_2(\sigma) = a_g(\sigma) - i b_g(\sigma). \quad (6) $$

This partition is very useful because the imaginary part, $b_g(\sigma)$, is directly related to the amount of energy dissipated by tides. It is therefore the main driver of the secular evolution of the system under the action of tides (see Sect. 2.3). For the Andrade model, we have from expression (5)

$$ b_g(\sigma) = k_f \frac{\mathcal{B}(\sigma) \sigma \tau_v}{\mathcal{A}(\sigma)^2 + \mathcal{B}(\sigma)^2} = (k_f - k_e) \frac{\mathcal{B}(\sigma) \sigma \tau}{\mathcal{A}(\sigma)^2 + \mathcal{B}(\sigma)^2}, \quad (7) $$

with

$$ \mathcal{A}(\sigma) = (\sigma\tau) \left[ 1 + |\sigma\tau|^{-\alpha} \left(\frac{\tau_e}{\tau}\right)^{1-\alpha} \cos\left(\frac{\alpha\pi}{2}\right) \Gamma(1+\alpha) \right], \quad (8) $$

and

$$ \mathcal{B}(\sigma) = 1 + |\sigma\tau|^{1-\alpha} \left(\frac{\tau_e}{\tau}\right)^{1-\alpha} \sin\left(\frac{\alpha\pi}{2}\right) \Gamma(1+\alpha), \quad (9) $$

where $\tau = \tau_e + \tau_v$ is the total relaxation time and $k_e = k_f \tau_e / \tau$ is the elastic Love number. For Earth and Venus, we have $k_e = 0.299$ and $k_e = 0.25$, respectively (Yoder 1995).

## 2.2. Thermal atmospheric tides

The differential absorption of the stellar heat by the planet's atmosphere gives rise to local variations of temperature and consequently to pressure gradients. The mass of the atmosphere is then permanently redistributed, adjusting for the equilibrium. As for gravitational tides, the redistribution of mass in the atmosphere gives rise to an atmospheric bulge that modifies the gravitational potential generated by the atmosphere in any point of the space, which is given by (e.g. Correia et al. 2003)

$$ V_a = -\hat{p}_2 \frac{3}{5\bar{\rho}} \left(\frac{R}{r}\right)^3 \left(\frac{a}{r_*}\right)^2 P_2(\hat{\mathbf{r}} \cdot \hat{\mathbf{r}}_*) \quad (10) $$

where $\bar{\rho}$ is the mean density of the planet and $\hat{p}_2$ is the second harmonic of the surface pressure variations, which depends on the frequency of the perturbation, $\sigma$, and on the composition and properties of the atmosphere.

In order to compute $\hat{p}_2(\sigma)$, we need to adopt some dynamical model for the atmosphere. For terrestrial planets in the habitable zone, the atmosphere can usually be treated as an ideal gas obeying the perfect gas law, together with the conservation of mass and the Navier-Stokes equations (e.g. Siebert 1961; Chapman & Lindzen 1970). In addition, assuming a slowly rotating convective atmosphere, the Coriolis acceleration can be neglected and the hydrodynamic equations greatly simplified

(e.g. Auclair-Desrotour et al. 2017a). Following Dobrovolskis & Ingersoll (1980), thermal tides are considered to be generated solely by the average heat absorbed at the ground, $J_0$, which gives for the second harmonic of the surface pressure variations (Auclair-Desrotour et al. 2017a, Eq. (166)):

$$ \hat{p}_2(\sigma) = -\frac{\kappa \rho_0 J_0}{\sigma_0 + i\sigma}, \quad (11) $$

with $\kappa = 1 - \gamma^{-1}$, where $\gamma$ is the adiabatic exponent value, $\rho_0$ is the mean surface density of the atmosphere, and $\sigma_0$ is the radiative cooling frequency of the atmosphere, which depends on its thermal capacity. As for the Love number (Eq. (6)), we can also decompose $\hat{p}_2$ in its real and imaginary parts as

$$ \hat{p}_2(\sigma) = a_a(\sigma) - i b_a(\sigma), \quad (12) $$

where the imaginary part characterises the viscous response of the atmosphere. The modulus gives the amplitude of the pressure variations and the argument gives the phase lag between the sub-stellar point and the maximal deformation of the atmosphere. We thus have from expression (11)

$$ b_a(\sigma) = -\frac{\kappa \rho_0 J_0 \sigma}{\sigma_0^2 + \sigma^2} = -\frac{\sigma/\sigma_0}{1 + (\sigma/\sigma_0)^2} \sqrt{\frac{10}{3\pi}} q_0, \quad (13) $$

where $q_0$ is the amplitude of the quadrupole term of the pressure variations at zero frequency (Leconte et al. 2015). This is the same functional form as in the model introduced by Ingersoll & Dobrovolskis (1978). The minus sign in the expression of $b_a(\sigma)$ causes the pressure variations to lead the star, a phenomenon well documented for Earth (e.g. Chapman & Lindzen 1970). In this study we do not consider the possible emerging effects of high-frequency resonances, such as the Lamb waves resonance or the effect of varying the thermal forcing profile (e.g. Farhat et al. 2024; Laskar et al. 2024).

## 2.3. Secular evolution

The additional disturbing potentials resulting from gravitational tides (Eq. (3)) and thermal tides (Eq. (10)) create a differential gravitational field around the planet given by

$$ g(\mathbf{r}) = -\nabla_r [V_g(\mathbf{r}) + V_a(\mathbf{r})]. \quad (14) $$

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

$$ \dot{G} = T = M_* \mathbf{r}_* \times g(\mathbf{r}_*), \quad (15) $$

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

$$ \dot{L} = -\dot{G} = -T. \quad (16) $$

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

$$ \dot{E} = M_* \dot{\mathbf{r}}_* \cdot g(\mathbf{r}_*). \quad (17) $$

In general, tidal effects slowly modify the spin and the orbit of the planet, on a timescale much longer than the orbital and precession periods of the system. We can then average the torque (Eq. (15)) and the power (Eq. (17)) over the mean anomaly and the argument of the pericentre, and obtain the equations of motion for the secular evolution of the system.

<footnote n="2">² The impact of $\tau_A \neq \tau_e$ has been studied in Bolmont et al. (2020).</footnote>
Following Correia & Valente (2022), we let (p, q, s) be a set of orthogonal vectors, such that

$$ \mathbf{p} = k \times \mathbf{s} \quad \text{and} \quad \mathbf{q} = s \times \mathbf{p} = k - x \mathbf{s}, \qquad (18) $$

with

$$ x = \cos \theta = k \cdot s, \qquad (19) $$

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. For the average torque, we therefore have

$$ \langle \mathbf{T} \rangle = (T_p^g + T_p^a) \mathbf{p} + (T_q^g + T_q^a) \mathbf{q} + (T_s^g + T_s^a) \mathbf{s}, \qquad (20) $$

and

$$ \langle \dot{\mathbf{E}} \rangle = n (T_E^g + T_E^a), \qquad (21) $$

where the superscript g refers to the contribution from gravitational tides (Sect. 2.3.1) and the superscript a refers to contribution from thermal atmospheric tides (Sect. 2.3.2). The Tp components in the torque (Eq. (20)) can be ignored because they change neither the norm of the angular momentum vectors nor the angle between them and therefore have no contribution to the secular spin and orbital evolution (see Sect. 2.4).

### 2.3.1. Gravitational tides

The complete expressions for the average torque and power due to gravitational tides, valid for any eccentricity, e, and obliquity, θ, are given by Correia & Valente (2022),

$$ T_q^g = \frac{3\mathcal{K}_g}{32} \sum_{k=-\infty}^{+\infty} \left\{ 3 b_g(-kn) (1-x^2) \left[ (X_k^{-3,2})^2 - (X_k^{-3,-2})^2 \right] \right. \\ \qquad \left. + 2 b_g(\omega - kn) \left[ (1+x)^2 (2-x) (X_k^{-3,2})^2 \right. \right. \\ \qquad \left. \left. - 4x^3 (X_k^{-3,0})^2 - (1-x)^2 (2+x) (X_k^{-3,-2})^2 \right] \right. \\ \qquad \left. + b_g(2\omega - kn) \left[ -4x(1-x^2)(X_k^{-3,0})^2 \right. \right. \\ \qquad \left. \left. + (1+x)^3 (X_k^{-3,2})^2 - (1-x)^3 (X_k^{-3,-2})^2 \right] \right\}, \qquad (22) $$

$$ T_s^g = \frac{3\mathcal{K}_g}{32} \sum_{k=-\infty}^{+\infty} \left\{ 2 b_g(\omega - kn) (1 - x^2) \left[ 4x^2 (X_k^{-3,0})^2 \right. \right. \\ \qquad \left. \left. + (1 + x)^2 (X_k^{-3,2})^2 + (1 - x)^2 (X_k^{-3,-2})^2 \right] \right. \\ \qquad \left. + b_g(2\omega - kn) \left[ 4(1 - x^2)^2 (X_k^{-3,0})^2 \right. \right. \\ \qquad \left. \left. + (1 + x)^4 (X_k^{-3,2})^2 + (1 - x)^4 (X_k^{-3,-2})^2 \right] \right\}, \qquad (23) $$

$$ T_E^g = \frac{\mathcal{K}_g}{64} \sum_{k=-\infty}^{+\infty} k \left\{ b_g(-kn) \left[ 4(1 - 3x^2)^2 (X_k^{-3,0})^2 \right. \right. \\ \qquad \left. \left. + 9(1 - x^2)^2 \left( (X_k^{-3,-2})^2 + (X_k^{-3,2})^2 \right) \right] \right. \\ \qquad \left. + 12 b_g(\omega - kn) (1 - x^2) \left[ 4x^2 (X_k^{-3,0})^2 \right. \right. \\ \qquad \left. \left. + (1 - x)^2 (X_k^{-3,-2})^2 + (1 + x)^2 (X_k^{-3,2})^2 \right] \right. \\ \qquad \left. + 3 b_g(2\omega - kn) \left[ 4(1 - x^2)^2 (X_k^{-3,0})^2 \right. \right. \\ \qquad \left. \left. + (1 - x)^4 (X_k^{-3,-2})^2 + (1 + x)^4 (X_k^{-3,2})^2 \right] \right\}, \qquad (24) $$

where

$$ \mathcal{K}_g = \frac{G M_*^2 R^5}{a^6}, \qquad (25) $$

and X_k^{-\ell,m}(e) are the Hansen coefficients, which depend only on the eccentricity<fnref n="3" /> (see Table 1).

### 2.3.2. Thermal atmospheric tides

The complete expressions for the average torque and power due to thermal tides are given by Valente & Correia (2023),

$$ T_q^a = \frac{3\mathcal{K}_a}{32} \sum_{k=-\infty}^{+\infty} \left\{ 3 b_a(-kn) (1-x^2) \left[ X_k^{-3,2} X_k^{-2,2} - X_k^{-3,-2} X_k^{-2,-2} \right] \right. \\ \qquad \left. + 2 b_a(\omega - kn) \left[ (1+x)^2 (2-x) X_k^{-3,2} X_k^{-2,2} \right. \right. \\ \qquad \left. \left. - 4x^3 X_k^{-3,0} X_k^{-2,0} - (1-x)^2 (2+x) X_k^{-3,-2} X_k^{-3,-2} \right] \right. \\ \qquad \left. + b_a (2\omega - kn) \left[ -4x(1-x^2) X_k^{-3,0} X_k^{-2,0} \right. \right. \\ \qquad \left. \left. + (1+x)^3 X_k^{-3,2} X_k^{-2,2} - (1-x)^3 X_k^{-3,-2} X_k^{-2,-2} \right] \right\}, \qquad (26) $$

$$ T_s^a = \frac{3\mathcal{K}_a}{32} \sum_{k=-\infty}^{+\infty} \left\{ 2 b_a(\omega - kn) (1-x^2) \left[ 4x^2 X_k^{-3,0} X_k^{-2,0} \right. \right. \\ \qquad \left. \left. + (1+x)^2 X_k^{-3,2} X_k^{-2,2} + (1-x)^2 X_k^{-3,-2} X_k^{-2,-2} \right] \right. \\ \qquad \left. + b_a(2\omega - kn) \left[ 4(1-x^2)^2 X_k^{-3,0} X_k^{-2,0} \right. \right. \\ \qquad \left. \left. + (1+x)^4 X_k^{-3,2} X_k^{-2,2} + (1-x)^4 X_k^{-3,-2} X_k^{-2,-2} \right] \right\}, \qquad (27) $$

$$ T_E^a = \frac{\mathcal{K}_a}{64} \sum_{k=-\infty}^{+\infty} k \left\{ b_a(-kn) \left[ 4(1 - 3x^2)^2 X_k^{-3,0} X_k^{-2,0} \right. \right. \\ \qquad \left. \left. + 9(1 - x^2)^2 (X_k^{-3,-2} X_k^{-2,-2} + X_k^{-3,2} X_k^{-2,2}) \right] \right. \\ \qquad \left. + 12 b_a(\omega - kn) (1 - x^2) \left[ 4x^2 X_k^{-3,0} X_k^{-2,0} \right. \right. \\ \qquad \left. \left. + (1 - x)^2 X_k^{-3,-2} X_k^{-2,-2} + (1 + x)^2 X_k^{-3,2} X_k^{-2,2} \right] \right. \\ \qquad \left. + 3 b_a(2\omega - kn) \left[ 4(1 - x^2)^2 X_k^{-3,0} X_k^{-2,0} \right. \right. \\ \qquad \left. \left. + (1 - x)^4 X_k^{-3,-2} X_k^{-2,-2} + (1 + x)^4 X_k^{-3,2} X_k^{-2,2} \right] \right\}, \qquad (28) $$

where

$$ \mathcal{K}_a = \frac{3M_* \left(\frac{R}{a}\right)^3}{5\bar{\rho}}. \qquad (29) $$

<footnote n="3">³ In this work we consider only terms with $|k| \le 22$, which gives an error of $e^{20} \sim 10^{-8}$ for the maximal adopted eccentricity ($e = 0.4$).</footnote>
<table><caption>**Table 1.** Hansen coefficients up to $e^6$.</caption><thead><tr><th>k</th><th>X<sub>k</sub><sup>-2,0</sup>(e)</th><th>X<sub>k</sub><sup>-2,2</sup>(e)</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>1223</sup>&frasl;<sub>320</sub>e<sup>6</sup></td><td></td><td><sup>3167</sup>&frasl;<sub>320</sub>e<sup>6</sup></td><td></td></tr><tr><td>-5</td><td><sup>1097</sup>&frasl;<sub>384</sub>e<sup>5</sup></td><td></td><td><sup>1773</sup>&frasl;<sub>256</sub>e<sup>5</sup></td><td></td></tr><tr><td>-4</td><td><sup>103</sup>&frasl;<sub>48</sub>e<sup>4</sup> - <sup>129</sup>&frasl;<sub>160</sub>e<sup>6</sup></td><td>-<sup>2</sup>&frasl;<sub>45</sub>e<sup>6</sup></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>13</sup>&frasl;<sub>8</sub>e<sup>3</sup> - <sup>25</sup>&frasl;<sub>128</sub>e<sup>5</sup></td><td>-<sup>27</sup>&frasl;<sub>640</sub>e<sup>5</sup></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>5</sup>&frasl;<sub>4</sub>e<sup>2</sup> + <sup>1</sup>&frasl;<sub>6</sub>e<sup>4</sup> + <sup>21</sup>&frasl;<sub>64</sub>e<sup>6</sup></td><td>-<sup>1</sup>&frasl;<sub>24</sub>e<sup>4</sup> - <sup>3</sup>&frasl;<sub>80</sub>e<sup>6</sup></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>e + <sup>3</sup>&frasl;<sub>8</sub>e<sup>3</sup> + <sup>65</sup>&frasl;<sub>192</sub>e<sup>5</sup></td><td>-<sup>1</sup>&frasl;<sub>24</sub>e<sup>3</sup> - <sup>13</sup>&frasl;<sub>384</sub>e<sup>5</sup></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>1</sup>&frasl;<sub>2</sub>e<sup>2</sup> + <sup>3</sup>&frasl;<sub>8</sub>e<sup>4</sup> + <sup>5</sup>&frasl;<sub>16</sub>e<sup>6</sup></td><td></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>e + <sup>3</sup>&frasl;<sub>8</sub>e<sup>3</sup> + <sup>65</sup>&frasl;<sub>192</sub>e<sup>5</sup></td><td>-e + <sup>3</sup>&frasl;<sub>8</sub>e<sup>3</sup> + <sup>7</sup>&frasl;<sub>192</sub>e<sup>5</sup></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>5</sup>&frasl;<sub>4</sub>e<sup>2</sup> + <sup>1</sup>&frasl;<sub>6</sub>e<sup>4</sup> + <sup>21</sup>&frasl;<sub>64</sub>e<sup>6</sup></td><td>1 - <sup>7</sup>&frasl;<sub>2</sub>e<sup>2</sup> + <sup>29</sup>&frasl;<sub>16</sub>e<sup>4</sup> - <sup>53</sup>&frasl;<sub>288</sub>e<sup>6</sup></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>13</sup>&frasl;<sub>8</sub>e<sup>3</sup> - <sup>25</sup>&frasl;<sub>128</sub>e<sup>5</sup></td><td>3e - <sup>69</sup>&frasl;<sub>8</sub>e<sup>3</sup> + <sup>369</sup>&frasl;<sub>64</sub>e<sup>5</sup></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>103</sup>&frasl;<sub>48</sub>e<sup>4</sup> - <sup>129</sup>&frasl;<sub>160</sub>e<sup>6</sup></td><td><sup>13</sup>&frasl;<sub>2</sub>e<sup>2</sup> - <sup>55</sup>&frasl;<sub>3</sub>e<sup>4</sup> + <sup>239</sup>&frasl;<sub>16</sub>e<sup>6</sup></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>1097</sup>&frasl;<sub>384</sub>e<sup>5</sup></td><td><sup>295</sup>&frasl;<sub>24</sub>e<sup>3</sup> - <sup>13745</sup>&frasl;<sub>384</sub>e<sup>5</sup></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>1223</sup>&frasl;<sub>320</sub>e<sup>6</sup></td><td><sup>345</sup>&frasl;<sub>16</sub>e<sup>4</sup> - <sup>10569</sup>&frasl;<sub>160</sub>e<sup>6</sup></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>69251</sup>&frasl;<sub>1920</sub>e<sup>5</sup></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>42037</sup>&frasl;<sub>720</sub>e<sup>6</sup></td><td></td><td><sup>73369</sup>&frasl;<sub>720</sub>e<sup>6</sup></td></tr></tbody></table>

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

## 2.4. Spin and orbital evolution

The set of Eqs. (20) and (21) allows us to track the secular evolution of the system using the angular momentum vectors and the orbital energy. Concerning the spin, the rotation rate is directly given from the rotational angular momentum (Eq. (2)),

$$ \omega = \frac{L \cdot s}{C} = \frac{\sqrt{L \cdot L}}{C}, \qquad (30) $$

while the angle between the orbital and equatorial planes (also known as obliquity) is obtained from both angular momentum vectors as (Eq. (19))

$$ \cos \theta = k \cdot s = \frac{G \cdot L}{\sqrt{(G \cdot G)(L \cdot L)}}. \qquad (31) $$

Concerning the orbit, the semi-major axis is directly given from the orbital energy (Eq. (1))

$$ a = -\frac{G M_{\star} m}{2E}, \qquad (32) $$

while the eccentricity is obtained from the orbital angular momentum (Eq. (1))

$$ e = \sqrt{1 - \frac{(G \cdot k)^2}{\beta^2 n^2 a^4}} = \sqrt{1 - \frac{G \cdot G}{\beta^2 n^2 a^4}}. \qquad (33) $$

## 3. Venus

Venus has a peculiar spin state among all terrestrial planets (e.g. Smith 1963; Carpenter 1964); it has a slow retrograde rotation, with a rotation period of 243.023 day ($\omega/n = 0.92462$) and an obliquity $\theta = 177.361^\circ$ (Margot et al. 2021). The present rotation

of Venus may represent an equilibrium between gravitational and thermal atmospheric tides (Gold & Soter 1969; Correia & Laskar 2001; Auclair-Desrotour et al. 2017a). It is thus a unique example in our Solar System of the tidal mechanisms described in Sect. 2, which can be used to benchmark our model.

### 3.1. Physical properties

The mass, radius, internal structure constant and Love numbers of Venus are well known (Table 2). The properties and composition of its atmosphere are also relatively well determined. Therefore, it is possible to use general circulation models to compute the amplitude of the torque induced by thermal tides (e.g. Leconte et al. 2013, 2015; Auclair-Desrotour et al. 2019; Wu et al. 2023). Indeed, in the case of Venus it has been shown that Eq. (13) complies with the thermal response of the atmosphere at various frequencies and it was estimated that $\sigma_0 = 23.8 \text{ yr}^{-1}$ and $q_0 = 2.01 \text{ mbar}$ (Table 2).

However, the rigidity and viscosity of Venus' mantle are not much constrained. Since Venus is similar to Earth in mass and size, we can adopt the Earth's parameters. Indeed, the rigidity and viscosity obtained for Venus using internal structure models with Andrade rheology are consistent with the values estimated for the Earth (Bolmont et al. 2020; Melini et al. 2022; Saliby et al. 2023). The rigidity of the Earth's upper mantle is $\mu \approx 80 \text{ GPa}$ (Karato & Wu 1993), but the whole mantle rigidity can be larger (see Table 2.1 in Sabadini et al. 2016). By volume-averaging the rigidity value we get $\mu \approx 176 \text{ GPa}$. The average effective viscosity is much more uncertain, $\eta \sim 10^{21} \text{ Pa s}$ (Karato & Wu 1993), which gives for the Maxwell relaxation time $\tau_e \sim 200 \text{ yr}$ (for more details see Bolmont et al. 2020). Actually, the viscosity may vary by several orders of magnitude, depending on the depth and temperature (e.g. Kirby & Kronenberg 1987). In the case of Earth, the surface post-glacial rebound due to the last glaciation about
<table><caption>**Table 2.** Physical and dynamical parameters adopted for the planets.</caption><thead><tr><th>Parameter</th><th>Venus</th><th>Earth</th><th>Kepler-1229 b</th></tr></thead><tbody><tr><td>M<sub>★</sub> (M<sub>☉</sub>)</td><td>1.00</td><td>1.00</td><td>0.54</td></tr><tr><td>m (M<sub>⊕</sub>)</td><td>0.82</td><td>1.00</td><td>2.74</td></tr><tr><td>R (R<sub>⊕</sub>)</td><td>0.95</td><td>1.00</td><td>1.40</td></tr><tr><td>a (au)</td><td>0.723</td><td>1.000</td><td>0.313</td></tr><tr><td>F<sub>★</sub> (W m<sup>-2</sup>)</td><td>2610</td><td>1366</td><td>606</td></tr><tr><td>p<sub>s</sub> (bar)</td><td>92.</td><td>1.</td><td>10.</td></tr><tr><td>q<sub>0</sub> (mbar)</td><td>2.01</td><td>11.8</td><td>26.9</td></tr><tr><td>σ<sub>0</sub> (yr<sup>-1</sup>)</td><td>23.8</td><td>145.</td><td>85.3</td></tr><tr><td>ξ</td><td>0.330</td><td>0.331</td><td>0.331</td></tr><tr><td>k<sub>e</sub></td><td>0.25</td><td>0.299</td><td>0.299</td></tr><tr><td>k<sub>f</sub></td><td>0.928</td><td>0.933</td><td>0.933</td></tr><tr><td>α</td><td>0.3</td><td>0.3</td><td>0.3</td></tr><tr><td>τ<sub>e</sub> (yr)</td><td>1468</td><td>1468</td><td>1468</td></tr></tbody></table>

**Notes.** The atmosphere parameters for Earth and Venus are from Leconte et al. (2015), and for Kepler-1229 b they are from Auclair-Desrotour et al. (2019). The inner structure parameters such as the internal structure constant, $\xi$; the elastic Love number, $k_e$; and the fluid Love number, $k_f$, are from Yoder (1995), and for Kepler-1229 b they are made equal to the Earth's values. The remaining parameters of the Kepler-1229 system are from Morton et al. (2016), where m is estimated from R, assuming an Earth-like density.

10⁴ yr ago is still ongoing, suggesting that the Earth's mantle relaxation time can be $\tau_e \sim 4000$ yr (Turcotte & Schubert 2002).

## 3.2. Dynamical constraints on Venus' mantle

The evolution of the rotation rate can be obtained from expressions (20) and (30) as

$$ \dot{\omega} = -\frac{1}{C} \langle T \rangle \cdot s = -\frac{T_s^g + T_s^a}{C}. \quad (34) $$

The current eccentricity of Venus ($e = 0.007$) and its spin state<fnref n="4" /> ($\omega_e/n = -0.92462$, $\theta_e = 2.639^{\circ}$), are well determined (Margot et al. 2021), where $\omega_e$ and $\theta_e$ are the current equilibrium rotation and obliquity, respectively. For simplicity, we can consider that Venus evolves in a circular orbit ($e \approx 0$) with zero obliquity ($\theta \approx 0$), and expression (34) simplifies as (Eqs. (23) and (27))

$$ \dot{\omega} = -\frac{3K_g}{2C} b_g (2\omega - 2n) - \frac{3K_a}{2C} b_a (2\omega - 2n). \quad (35) $$

Assuming that the present rotation of Venus is in an equilibrium state<fnref n="5" />, we must have $\dot{\omega} = 0$ for $\omega = \omega_e$, and thus

$$ K_g b_g (2\omega_e - 2n) = -K_a b_a (2\omega_e - 2n). \quad (36) $$

The constant parameters $K_g$ (Eq. (25)) and $K_a$ (Eq. (29)), along with the surface pressure variation function, $b_a(2\omega_e - 2n)$ (Eq. (13)), are well determined for Venus (Table 2). The major

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Maxwell relaxation time $\tau_e$ as a function of the $\alpha$ parameter for Venus (Eq. (36)) obtained for the present rotation state.</figcaption></figure>

source of uncertainty is related to the Love number function, $b_g(2\omega_e - 2n)$, which depends on the mantle rheology. For the Andrade model (Sect. 2.1), this uncertainty is transferred to the determination of the $\alpha$ and $\tau_e$ parameters (Eq. (7)).

For materials deforming according to an Andrade rheology with linear work hardening, we have $\alpha = 1/3$ (Louchet & Duval 2009). However, for rocky planets it is commonly accepted that $0.2 < \alpha < 0.4$ (Castillo-Rogez et al. 2011). Indeed, for Earth, Tobie et al. (2019) have shown that $\alpha \in [0.2, 0.28]$ is able to correctly reproduce the observed tidal dissipation rate. Melini et al. (2022) have studied the tidal deformation of Venus with the ALMA³ code (plAnetary Love nuMbers cAlculator), and conclude that $\alpha = 0.3$ provides the best agreement with tidal constraints on the interior of Venus determined by Dumoulin et al. (2017).

In the case of Venus, we can use the dynamical constraint provided by the present spin state to find a link between $\alpha$ and $\tau_e$. In Fig. 1, we show the solution of Eq. (36) for these two parameters. We observe that $\tau_e$ can vary by several orders of magnitude for the expected range of $\alpha \in [0.2, 0.4]$. For $\alpha = 0.2$ it is about 60 kyr, while for $\alpha = 0.4$ it is only about 200 yr. Following Melini et al. (2022), we adopt $\alpha = 0.3$, which gives for the Maxwell relaxation time, $\tau_e = 1468$ yr. These values are consistent with those of the Earth (Sect. 3.1), so we adopted them for this class of planets in the remainder of this work (see Table 2).

## 3.3. Equilibrium states

Using the tidal models for Venus' mantle and atmosphere that better match the current observations (Sects. 3.1 and 3.2), we can quickly investigate all the possible equilibria for its spin.

In Fig. 2 a, we show the evolution of the rotation rate for zero obliquity (Eq. (35)). We observe that there are five equilibrium points ($\dot{\omega} = 0$), though only three correspond to the stable equilibria (those with $\partial\dot{\omega}/\partial\omega < 0$). One of these points corresponds to $\omega = n$ (synchronous rotation), while the remaining two correspond to the present state and its symmetric, $\omega - n = \pm\omega_s$, with $\omega_s/n = 1.92$ (Correia & Laskar 2001). The current observed state, $\omega_e = n - \omega_s$, is marked with a dot. The unstable equilibria are located at $\omega - n = \pm\omega_c$, with $\omega_c/n = 1.31$. These unstable points are also important to understand the behaviour of the rotation rate. For $\omega/n \in [-0.31, 2.31]$, the rotation will evolve into the synchronous state, while outside this interval it will move to the closest asynchronous stable state, $\omega = n \pm \omega_s$. These results are in agreement with previous studies that used

<footnote n="4">⁴ The current spin state of Venus can be described either by the pair $(\omega, \theta)$ or by the pair $(-\omega, \pi - \theta)$. These two points do not correspond to the same physical state, but they are equivalent from a dynamical point of view (Correia & Laskar 2001).</footnote>

<footnote n="5">⁵ Because Venus is close to synchronous rotation, it is often assumed that it already reached an equilibrium state, although this hypothesis is yet to be confirmed (e.g. Margot et al. 2021; Revol et al. 2023).</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Variation of $d\omega/dt$ upon $\omega/n$ for Venus with $e = 0$ (Eq. (34)). We show the evolution for different obliquity values: $\theta = 0^\circ$ (top), $\theta = 90^\circ$ (middle), and $\theta = 180^\circ$ (bottom plot). The black dots represent the presently observed rotation state.</figcaption></figure>

different tidal models (Correia & Laskar 2001, 2003; Leconte et al. 2015; Auclair-Desrotour et al. 2017b).

Because our model is valid for any obliquity value (Sect. 2), we can study the evolution of the rotation rate at any obliquity (Eq. (34)). In Fig. 2 c, we show the evolution for $\theta = 180^\circ$. As expected, this case is similar to the case with $\theta = 0^\circ$, since the pair ($\omega, 180^\circ$) behaves like the pair ($-\omega, 0^\circ$) (Correia & Laskar 2001). However, assuming that Venus rotation rate brakes from initial fast rotations, we always come from the right-hand side of the figure. As a consequence, for $\theta = 0^\circ$, Venus is expected to evolve into the state $\omega = n + \omega_s$, while for $\theta = 180^\circ$, it is expected to end up in the state $\omega = n - \omega_s$ (present state).

The previous two obliquity values ($\theta = 0^\circ$ and $180^\circ$) are the most likely outcomes of gravitational and thermal tides (Correia & Laskar 2001, 2003; Correia et al. 2003). However, we cannot rule out that other possibilities for the obliquity exist, at least temporarily. Therefore, in Fig. 2 b, we also show the evolution for $\theta = 90^\circ$. Interestingly, we observe that the asynchronous thermal states are no longer possible. On the other hand, the spin evolution is dominated by spin-orbit resonances that were absent before ($\omega/n = -2, -1, 0, +1, +2$). These resonant states result from the geometry of the tidal torques (Eqs. (7) and (13)), which experience a significant increment at these specific rotation rates (Eqs. (23) and (27)). As a result, the spin will settle into one of these states and stay there until the obliquity evolves into more modest values. Although $\theta = 90^\circ$ may not correspond to a final possibility for the spin of Venus, the crossing of the spin-orbit resonances can modify its final evolution (see next section).

### 3.4. Complete spin evolution

The rotation rate and the obliquity evolution cannot be dissociated (Eqs. (30) and (31)). Therefore, the exact evolution of the

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Secular evolution of the spin of Venus in the plane $(\omega/n, \theta)$ for different initial obliquity values. This plot was obtained with an initial $\omega/n = 4.5$ and by integrating Eq. (16), where the Hansen coefficients $X_k^{l,m}(e)$ are truncated at $|k| \le 22$. The stable equilibrium points are marked with a black dot.</figcaption></figure>

spin can only be obtained by integrating Eq. (16). This allows us to check how Venus may have evolved into the different equilibrium states (Sect. 3.3) and also the existence of some eventual nonzero obliquity equilibria.

The orbital evolution (Eqs. (15) and (17)) is negligible in the case of Venus because $||L|| \ll ||G||$. Actually, the eccentricity of Venus undergoes much stronger variations due to planetary perturbations than due to tidal effects (Correia & Laskar 2003). The present eccentricity of Venus is about 0.007 and its mean value is approximately 0.022 (Laskar 2008). For simplicity, we thus consider that Venus evolves in a circular orbit ($e \approx 0$) and that G and E are constant (Eq. (1)).

In Fig. 3, we show some trajectories for the tidal evolution of the spin in the plane $(\omega/n, \theta)$. We start the integrations with $\omega/n = 4.5$, and different initial obliquity values in the range $[0^\circ, 180^\circ]$. Tidal effects initially decrease the rotation rate (Eq. (34)), so the evolutionary paths must be followed for decreasing values of $\omega/n$, until they reach an equilibrium position (marked with a black dot).

For the initial low obliquity ($\theta = 0^\circ$), the obliquity does not change much and the rotation rate reaches the asynchronous rotation with $\omega/n = 2.92$, corresponding to the prograde final state ($\omega = n + \omega_s, \theta = 0^\circ$). For initial high obliquity ($\theta \ge 160^\circ$), the obliquity evolves into $180^\circ$ and then the rotation rate reaches the asynchronous rotation with $\omega/n = 0.92$, which is equivalent to the present retrograde final state ($\omega = n - \omega_s, \theta = 0^\circ$).

For initial obliquities ranging from $20^\circ$ to $140^\circ$, the rotation rate can be locked in the 2/1 spin-orbit (see also Fig. 2 b). At that point, the obliquity decreases until it reaches a much smaller value. At the critical value $\theta \approx 3^\circ$, the 2/1 resonance becomes unstable and the rotation rate is allowed to decrease again. It then moves to the synchronous rotation, corresponding to the final state ($\omega = n, \theta = 0^\circ$). This behaviour can be perfectly understood by looking at Fig. 2 a. For near zero obliquities, rotation rates
larger than the critical value $\omega/n = 2.31$ move into the asynchronous final state $\omega = n + \omega_s$, while rotation rates smaller than these critical value move into the synchronous state.

Our results are in perfect agreement with those obtained by Revol et al. (2023), who studied the spin evolution of Venus-like planets using different values for the $\alpha$ and $\tau_e$ parameters. However, our results do not explain so well the current observations, since the present state can only be reached by a small set of initially extremely high obliquities ($\theta \ge 160^{\circ}$). This is likely because in the present study we are not including all the effects that contribute to modify the spin evolution of Venus. Indeed, we do not take into account either the effect of core-mantle friction, which tends to align the equator of the planet with the orbital plane (Correia et al. 2003), or the effect of planetary perturbations, which can introduce chaotic variations on the obliquity (Correia & Laskar 2001, 2003).

## 4. Earth-like planets

A configuration similar to the Earth around the Sun is ideal when searching for habitable worlds. However, low-mass stars are more frequent than sun-like stars in the Milky Way (e.g. Bochanski et al. 2010). The occurrence rate of planets around low-mass stars is also higher because they are easier to detect (e.g. Winn 2018). Moreover, formation studies show that temperate Earth-sized planets may be more frequent around M and K-dwarf stars (Burn et al. 2021). The fact that Earth-like planets can be relatively common around low-mass stars offers a great opportunity to explore the possibility of having life.

The mean surface pressure on Earth is only $p_s = 1$ bar, which contrasts with that of Venus, $p_s = 92$ bar (Table 2). However, the amplitude of the thermal tide for an Earth-like atmosphere is about an order of magnitude larger, because on Venus the sunlight is almost completely scattered or absorbed before it reaches the surface (Leconte et al. 2015). As a result, Earth-like planets around low-mass stars can develop non-synchronous rotation similar to that of Venus (Correia et al. 2008; Cunha et al. 2015; Leconte et al. 2015).

In this section, we revisit the equilibria for the spin of Earth-like planets around low-mass stars. We benefited from more accurate models for both gravitational and thermal atmospheric tides (Sect. 2), which allows us to perform a more comprehensive study of all the possibilities for a wide range of stellar masses, semi-major axes, and eccentricities. We first explain how we obtain the stellar mass as a function of the semi-major axis that preserves the average flux received by the planet (Sect. 4.2). We then determine the number of possible equilibrium states for a given semi-major axis (Sect. 4.3) and show how the position of the thermal asynchronous non-resonant equilibria change with the semi-major axis and with the eccentricity (Sect. 4.4). Finally, we generalize these results for any planet in the habitable zone, that is, with an arbitrarily average stellar flux (Sect. 4.5).

### 4.1. Physical properties

We consider an Earth-twin planet around a low-mass star, that is, we assume the exact same mass, radius, atmosphere composition and average annual stellar flux as the Earth's (Table 2). It is then possible to run general circulation models to compute the amplitude of the torque induced by thermal tides, which was done in Leconte et al. (2015). For an average flux $F_\star = 1366$ W m$^{-2}$,

they obtain<fnref n="6" /> for the radiative frequency of the atmosphere, $\sigma_0 = 145$ yr$^{-1}$, and for the amplitude of the atmospheric quadrupole, $q_0 = 11.8$ mbar. For the gravitational tides we adopt the Andrade model with the parameters $\alpha = 0.3$ and $\tau_e = 1468$ yr (Sect. 3.2).

### 4.2. Stellar mass

As we fix the average incoming stellar flux, we can determine the stellar luminosity from the orbital parameters of the planet through the following relation (e.g. Laskar et al. 1993)

$$F_{\star} = \left\langle \frac{L_{\star}}{4\pi r_{\star}^{2}} \right\rangle = \frac{L_{\star}}{4\pi a^{2} \sqrt{1 - e^{2}}}. \qquad (37)$$

We then use the mass-luminosity relation to obtain the stellar mass. For the low-mass stars under investigation, this relation is obtained for $0.20 < M_\star/M_\odot < 0.85$ as (Cuntz & Wang 2018),

$$\frac{L_{\star}}{L_{\odot}} \approx \left( \frac{M_{\star}}{M_{\odot}} \right)^{f(M_{\star}/M_{\odot})}, \qquad (38)$$

with

$$f(x) = -141.7x^4 + 232.4x^3 - 129.1x^2 + 33.29x + 0.215, \qquad (39)$$

where $L_\odot$ and $M_\odot$ are the Sun luminosity and mass, respectively. By solving Eqs. (37) and (38), we get the stellar mass as a function of the semi-major axis and eccentricity for a constant flux.

### 4.3. Equilibrium final states

In the absence of thermal atmospheric tides, the torque resulting from gravitational tides (Sect. 2.3.1) drives the spin into synchronous rotation, $\omega/n = 1$, and zero obliquity, $\theta = 0^{\circ}$ (e.g. Hut 1980). In the presence of thermal tides, an additional torque on the spin emerges (Sect. 2.3.2) that can counter the braking effect from gravitational tides and give rise to additional equilibrium possibilities. In order to investigate the multitude of scenarios, in Fig. 4, we show the evolution of the rotation rate with zero obliquity for different values of the semi-major axis and eccentricity (Eq. (34)),

$$\dot{\omega} = \frac{3}{2C} \sum_{k=-\infty}^{+\infty} \left[ \mathcal{K}_g b_g (2\omega - kn) (X_k^{-3,2})^2 + \mathcal{K}_a b_a (2\omega - kn) X_k^{-3,2} X_k^{-2,2} \right]. \qquad (40)$$

Equilibrium points are observed for $\dot{\omega} = 0$, and they can be stable for negative slopes ($\partial\dot{\omega}/\partial\omega < 0$).

For circular orbits (Fig. 4, red curves), there is a critical semi-major axis ($a_c \approx 0.32$ au) below which the spin behaves as in absence of thermal tides, that is, there is only one equilibrium possibility, corresponding to the synchronous rotation (Fig. 4a). For an Earth-like planet, this critical threshold corresponds to a stellar mass $M_\star/M_\odot \approx 0.68$. However, for larger semi-major axes (hence, for larger stellar masses), we observe that there are five equilibrium points (Fig. 4c), as in the case of Venus (Sect. 3.3).

<footnote n="6">⁶ The response of the atmosphere in Leconte et al. (2015) is obtained using general circulation models for a planet around a Solar-type star. We note that, eventually, the same average stellar flux generated by a different stellar type star may not be absorbed exactly in the same way.</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Variation of $d\omega/dt$ upon $\omega/n$ for an Earth-like planet with $\theta = 0^\circ$ (Eq. (40)). We show the evolution for different semi-major axes: $a/a_c = 0.63$ ($M_*/M_\odot \approx 0.52$; top), $a = a_c = 0.32$ au ($M_*/M_\odot \approx 0.68$; middle), and $a/a_c = 1.10$ ($M_*/M_\odot \approx 0.70$; bottom plot). We also show the evolution for three different values of the eccentricity: $e = 0.0$ (red line), $e = 0.2$ (orange line), and $e = 0.4$ (blue line).</figcaption></figure>

Again, three of these points can be stable, one for $\omega = n$ (syn-
chronous rotation) and the other two for $\omega = n \pm \omega_s$. Here,
$\omega_s = |\omega - n|$ is the synodic rotation frequency (that is, the fre-
quency between two passages of the star in the local meridian),
which depends on the balance between thermal and gravitational
tides. We thus confirm that the rotation of Earth-like planets in
circular orbits around low-mass stars can also be locked in asyn-
chronous states (Leconte et al. 2015). The emergence of these
new equilibria is favoured for planets further away from the star,
because the relative strength between thermal and gravitational
tides increases for larger semi-major axis (Correia et al. 2008).
Indeed, from expressions (25) and (29), we see that $K_d/K_g \propto a^3$.

Gravitational tides tend to circularise the orbits, but usually
on a much longer timescale than they take to modify the spin
(e.g. Valente & Correia 2022). Moreover, if the system has sev-
eral planets, their mutual gravitational interactions permanently
excite the eccentricities (e.g. Laskar 1994, 2008). Therefore, we
expect that Earth-like planets can be observed with some ecen-
tricity, provided that the system remains stable. For eccentric
orbits (Fig. 4, orange and blue curves), the number of possible
equilibria for the spin increases and the figure becomes more
complex. In this case, gravitational tides alone can already drive
the rotation rate into asynchronous spin-orbit resonances (e.g.
Makarov & Efroimsky 2013; Correia et al. 2014). These new
states centred at $\omega/n = k/2$ (with $k \in \mathbb{Z}$) arise from the contri-
bution of the many Hansen coefficients (Eq. (40) and Table 1).
Indeed, higher-order spin-orbit resonances appear as we increase
the eccentricity. When we take into account thermal tides, addi-
tional nonzero Hansen coefficients are present at the exact same
frequencies (Eq. (40)), but do not change much the balance of
torques. Gravitational tides thus control the spin evolution near

the spin-orbit resonances, which remain a stable possible out-
come for the spin evolution. Nevertheless, as in the case with
zero eccentricity, thermal tides can still introduce an additional
asynchronous and non-resonant state at both ends of the spin-
orbit resonance region, although these equilibria become more
difficult to attain (the critical semi-major axis $a_c$ increases with
the eccentricity). In the example shown in Fig. 4c (with $a =
0.35$ au), for $e = 0.2$, the thermal asynchronous equilibria are
still present, but for $e = 0.4$, they are no longer possible.

4.4. Asynchronous non-resonant thermal states

During the formation stages, planets are expected to acquire fast
rotation rates (e.g. Kokubo & Ida 2007), such that we initially
have $|\omega| \gg n$. Then, for zero obliquity, as the rotation approaches
the equilibria region ($\omega \sim n$), it first encounters and locks in one
of the asynchronous non-resonant equilibria introduced by ther-
mal tides (provided that these states exists for a given semi-major
axis and eccentricity). Since asynchronous thermal states are of
major interest for habitability, it is important to better understand
how they depend on the orbital parameters.

In Fig. 5, we first determine the position of the asynchronous
thermal states as a function of the semi-major axis, for some
fixed values of the eccentricity ($e = 0, e = 0.2, e = 0.3$, and
$e = 0.4$). The different equilibria are obtained numerically by
solving $\dot{\omega} = 0$ (Eq. (40)). We show the corresponding relative
synodic frequency, $\omega_s/n = |\omega/n - 1|$, but also the synodic rota-
tion period, $P_s = 2\pi/\omega_s$ (in days). The solid lines follow the
state with $\omega/n > 1$ (super-synchronous), while the dashed lines
map the state with $\omega/n < 1$ (sub-synchronous). The synodic
frequency is the same for zero eccentricity (Correia & Laskar
2001), but it can be slightly different for eccentric orbits, in par-
ticular close to the critical semi-major axis. We recall that in
our analysis we fix the stellar flux at the same value as that of
the Earth's. Therefore, by modifying the semi-major axis, we
also change the stellar mass according to expressions (37) and
(38). The mass range $0.20 < M_*/M_\odot < 0.85$ corresponds to a
semi-major axis range $0.07 < a < 0.72$ au.

In Fig. 5, we observe that the relative equilibrium rotation
rate, $\omega/n$, increases with the semi-major axis (and hence with the
stellar mass). Conversely, the synodic rotation period decreases
with the semi-major axis, that is, the planet rotates faster. Inter-
estingly, we find that the equilibrium can be close to $P_s \approx 1$ day
(the current synodic period of the Earth), for some semi-major
axes. More precisely, for $a \approx 0.62$ au when $e = 0$ (circular orbit),
for $a \approx 0.64$ au when $e = 0.2$ (small eccentricity), for $a \approx 0.66$ au
when $e = 0.3$, and for $a \approx 0.70$ au when $e = 0.4$ (moderate
eccentricity). These semi-major axes provide the same stellar
flux of the Earth for a star with $M_* \approx 0.82 M_\odot$. We hence con-
clude that Earth-like planets around low-mass stars can maintain
a rotation period similar to that of the current Earth despite being
closer to their parent stars.

In Fig. 5, we can also clearly see that there is a critical
semi-major axis below which the asynchronous thermal states
no longer exist. Planets with semi-major axes below the critical
values thus evolve into some spin-orbit resonance (Fig. 4).
We further confirm that the critical semi-major axes depend
on the eccentricity. As the eccentricity increases, the criti-
cal value of the semi-major axis also increases (check also
Fig. 4c). Moreover, for eccentric orbits, we can spot a small
difference in the critical semi-major axes between the sub- and
super-synchronous equilibrium states. The equilibrium rotations
between these two states also differ more significantly close to
the critical semi-major axes. These subtle differences become
more visible for the case with higher eccentricity ($e = 0.4$).
<figure><img src="image_6.png" /></figure>

**Fig. 5.** Relative synodic frequency, $\omega_s/n = |\omega/n - 1|$ (top), and respective synodic rotation period, $P_s = 2\pi/\omega_s$ (bottom), of the stable asynchronous thermal equilibrium states as a function of the semi-major axis for different values of the eccentricity, $e = 0, 0.2, 0.3, 0.4$. The solid lines give the state with $\omega/n > 1$ (super-synchronous), while the dashed lines give the state with $\omega/n < 1$ (sub-synchronous). These plots are shown on a log-scale in the y-axis.

Since the eccentricity introduces some differences in the final equilibria, in Fig. 6, we determine the position of the asynchronous thermal states as a function of the eccentricity, for some fixed values of the semi-major axis from $a = 0.35$ au to 0.6 au. All these semi-major axis values are larger than the critical semi-major axis for a circular orbit. Because we fix the stellar flux, the stellar masses are more or less<fnref n="7" /> constant for each semi-major axis, going from $M_*/M_\odot \approx 0.7$ to 0.8, respectively. We show again the relative synodic frequency, $\omega_s/n = |\omega/n - 1|$, together with the corresponding synodic rotation period,

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

**Fig. 6.** Relative synodic frequency, $\omega_s/n = |\omega/n - 1|$ (top), and respective synodic rotation period, $P_s = 2\pi/\omega_s$ (bottom), of the stable asynchronous thermal equilibrium states as a function of the eccentricity for different values of the semi-major axis, $a = 0.35, 0.4, 0.5, 0.6$ au, corresponding to the stellar masses $M_*/M_\odot \approx 0.70, 0.73, 0.77, 0.81$, respectively. The solid lines give the state with $\omega/n > 1$ (super-synchronous), while the dashed lines give the state with $\omega/n < 1$ (sub-synchronous). These plots are shown on a log-scale in the y-axis.

$P_s = 2\pi/\omega_s$ (in days). The solid lines show the state with $\omega/n > 1$ (super-synchronous), while the dashed lines give the state with $\omega/n < 1$ (sub-synchronous). We confirm that the synodic frequency is the same in the two states for zero eccentricity, but the difference increases with the eccentricity, in particular close to the critical semi-major axis. We also observe that the effect of the eccentricity is more important for smaller semi-major axis. For $a = 0.35$ au, the critical semi-major axis is attained for $e \approx 0.33$, which explains why no asynchronous thermal equilibria is observed for $e = 0.4$ in Fig. 4c.

### 4.5. Habitable zone

In the previous section, we have seen that for Earth-like planets, asynchronous thermal states can only exist for semi-major

<footnote n="7">⁷ Actually, the stellar flux also changes with the eccentricity (Eq. (37)), and hence the stellar mass. However, the impact of the eccentricity is not very important, because for the maximum value chosen ($e = 0.4$), we still have $\sqrt{1 - e^2} \approx 0.92$.</footnote>
axes larger than some critical value, $a_c \approx 0.35$ au, which correspond to stars with $M_*/M_\odot \ge 0.7$ (K-dwarf stars). However, this analysis was restricted to an exact Earth-twin planet, that is, a planet with the same mass, radius, atmosphere composition, and average annual stellar flux as the Earth's (Table 2).

The existence of life as we know it, based on liquid water, is likely not so restrictive. In principle, habitable worlds can exist for a board range of stellar fluxes, $450 < F_* < 1508$ Wm$^{-2}$, known as habitable zone (e.g. Kasting et al. 1993; Leconte et al. 2015). The amplitude of thermal tides depends on the stellar flux (Eq. (11)), so the position of the asynchronous thermal states will also change. For an Earth-like planet ($N_2$ atmosphere and $p_s = 1$ bar), Leconte et al. (2015) also compute the thermal response for the stellar flux $F_* = 450$ Wm$^{-2}$ using general circulation models. They obtain for the radiative frequency of the atmosphere, $\sigma_0 = 74.5$ yr$^{-1}$, and for the amplitude of the atmospheric quadrupole, $q_0 = 8.90$ mbar. It is then possible to determine the position of the different equilibria for the spin for this specific stellar flux. We can also determine the critical semi-major axis. For a circular orbit we get $a_c = 0.34$ au, which corresponds to a stellar mass $M_*/M_\odot = 0.51$. This critical semi-major axis is slightly larger than the one obtained for Earth's stellar flux (Fig. 4 b), but the stellar mass is much smaller and already in the M-dwarf regime. Since the critical semi-major axes appears to not be very sensitive to the stellar flux, from the two points provided by Leconte et al. (2015) at $F_* = 450$ Wm$^{-2}$ and $1366$ Wm$^{-2}$, one can extrapolate using a straight line that, for Earth-like planets, we have

$$a_c \approx (0.4 - 0.12 M_*/M_\odot) \text{ au.} \quad (41)$$

In order to get a more comprehensive view of the possible behaviour of Earth-like planets at different stellar fluxes, in Fig. 7, we show the critical semi-major axes (Eq. (41)) as a function of the stellar mass and semi-major axis for circular orbits (blue line). Only planets on the right-hand side of this line can be found in an asynchronous thermal state, while those on the left-hand side will likely develop synchronous rotations (at least for zero eccentricity). In Fig. 7, we also show the limits of the habitable zone (black lines), to better guide us to the most interesting areas. The top line corresponds to $F_* = 1508$ Wm$^{-2}$, while the bottom line corresponds to $F_* = 450$ Wm$^{-2}$. We observe that for $M_*/M_\odot \ge 0.5$ and $a \ge 0.35$ au, it is possible to find a wide region where planets can be simultaneously asynchronous and inside the habitable zone<fnref n="8" />.

Habitable worlds can also exist for a large variety of planetary properties (e.g. Kopparapu et al. 2013, 2020), that is, they are not restricted to the Earth's mass, radius, and atmosphere. Actually, because of the current detection limitations, most of the low-mass (rocky) planets known to date are in the super-Earth regime (e.g. Schlichting 2018). For comparison, in Fig. 7, we plot all the known planets with $R/R_⊕ < 2.2$ (which is equivalent to $m/M_⊕ < 10$, assuming a mean density identical to that of the Earth's). For planets with $R/R_⊕ < 1.7$ (that is, with $m/M_⊕ < 5$, marked with a black dot), the mass is not large enough to accrete a voluminous gaseous envelope, so they are likely composed of a rocky mantle and a thin atmospheric layer, as observed for Earth

<figure><img src="image_8.png" /><figcaption>Fig. 7. Critical semi-major axes for Earth-like planets around low-mass stars. We show the results for an atmosphere with $p_s = 1$ bar (in blue) and $p_s = 10$ bar (in red). Planets on the right-hand side of the critical lines are expected to present asynchronous rotation. We also show the limits of the habitable zone (black curves) for stellar fluxes within $450 < F_* < 1508$ Wm$^{-2}$. The dots correspond to the already known planets (in black for $R/R_⊕ < 1.7$, in green for $1.7 < R/R_⊕ < 2.0$, and in orange for $2.0 < R/R_⊕ < 2.2$).</figcaption></figure>

and Venus (e.g. Demory et al. 2016; Komacek & Abbot 2019; Wordsworth & Kreidberg 2022).

For rocky planets in the super-Earth regime, we nevertheless expect a more dense atmosphere than that of the Earth. Using general circulation models, Auclair-Desrotour et al. (2019) derived some scaling laws for $q_0$ and $\sigma_0$ as a function of the semi-major axis, for a dry terrestrial planet with a homogeneous $N_2$ atmosphere and $p_s = 10$ bar around a Sun-like star<fnref n="9" />. These laws can be easily converted in terms of stellar flux as (Eq. (37))

$$\log_{10} a = 1.568 - \frac{1}{2} \log_{10} F_{*}, \quad (42)$$

which gives

$$\log_{10} q_0 = 0.345 \log_{10} F_* + 2.469 \quad (43)$$

and

$$\log_{10} \sigma_0 = 0.430 \log_{10} F_* + 0.7343, \quad (44)$$

where $a$ is in au, $F_*$ is given in Wm$^{-2}$, $q_0$ in Pa, and $\sigma_0$ in yr$^{-1}$. Therefore, it is possible to know the thermal tidal response of the atmosphere for any planet of this kind, and hence the position of the asynchronous thermal states.

In Fig. 7, we additionally show the critical semi-major axes for an Earth-like planet with $p_s = 10$ bar (red line). Interestingly, we observe that the critical semi-major axis and the stellar mass required to prevent synchronous rotation are much smaller than

<footnote n="8">⁸ Our Fig. 7 is similar to Fig. 3 in Leconte et al. (2015), but it has some important differences for the critical semi-major axes. They used the constant-$Q$ model for gravitational tides, while we adopted the Andrade model (Sect. 2.1). More importantly, in Eq. (41) we compute the mass of the star that is able to provide a given stellar flux $F_*(a_c, M_*)$ according to expressions (37) and (38).</footnote>

<footnote n="9">⁹ From Table 2 in Auclair-Desrotour et al. (2019), the scaling laws are $\log_{10} q_m = -0.69 \log_{10} a + 3.25$, and $\log_{10} \tau_m = 0.86 \log_{10} a + 0.48$, where $a$ is given in au, $q_m$ in Pa, and $\tau_m$ in day, $q_0 = 2q_m$ and $\sigma_0 = 1/\tau_m$.</footnote>
for a planet with $p_s = 1$ bar (blue line). Many of the already known planets are in the right-hand side of this critical line, so we expect them to present asynchronous rotation. Some of these planets are also inside the habitable zone, which can foster friendly conditions for life to develop.

## 5. Application to Kepler-1229 b

In order to study the complete long-term tidal evolution of a given Earth-like system over time, we need to attribute specific values to all physical and dynamical parameters. One good candidate for habitability studies is Kepler-1229 b (Morton et al. 2016), a super-Earth with a radius of $R = 1.40 R_⊕$ in a 87 day orbit around a star with mass $M_★ = 0.54 M_⊙$ ($a = 0.3125$ au), which is in the middle of the habitable zone (Fig. 7). The eccentricity is not constrained because so far this system was solely observed using the transits technique.

Despite Kepler-1229 b being slightly larger than Earth, we adopted here a similar average density, inner structure constant, Love numbers, relaxation time, and Andrade rheology. As a result, we obtained a larger mass, $m \approx 2.74 M_⊕$, which is able to accrete more gas. Therefore, we assumed a denser $N_2$ atmosphere for Kepler-1229 b, with $p_s = 10$ bar. For a circular orbit<fnref n="10" />, we get from Eq. (37) that $F_★ \approx 606$ Wm⁻², which in turn allowed us to obtain for the amplitude of the atmospheric quadrupole, $q_0 = 26.9$ mbar (Eq. (43)), and for the radiative frequency of the atmosphere, $\sigma_0 = 85.3$ yr⁻¹ (Eq. (44)). In Table 2, we summarise the physical and dynamical parameters that we used.

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Variation of $d\omega/dt$ upon $\omega/n$ for Kepler-1229 b. We show the evolution for different obliquity values: $\theta = 0°$ (top), $\theta = 45°$ (middle), and $\theta = 90°$ (bottom plot). We also show the evolution for three different values of the eccentricity: $e = 0.0$ (red line), $e = 0.2$ (orange line), and $e = 0.4$ (blue line).</figcaption></figure>

### 5.1. Equilibrium final states

We first investigate the possible equilibrium final states for the given parameters of Kepler-1229 b. In Fig. 8a, we show the evolution of the rotation rate for zero obliquity (Eq. (40)) and different values of the eccentricity ($e = 0, e = 0.2, \text{and } e = 0.4$). As for the case of Venus (Sect. 3.3) and for an Earth-like planet with $a > a_c$ (Sect. 4.3), we observe that for zero eccentricity there are three stable equilibrium points, one for $\omega = n$ (synchronous rotation), and another two for $\omega = n \pm \omega_s$ (asynchronous thermal states). More precisely, we find $\omega_s/n \approx 4.1$ (corresponding to $P_s \approx 21$ day), which gives one prograde state at $\omega/n \approx 5.1$, and one retrograde state at $\omega/n \approx -3.1$. As we increase the eccentricity, additional equilibria appear at spin-orbit resonances $\omega/n = k/2$ (with $k \in \mathbb{Z}$), while the asynchronous thermal states tend to disappear. They still exist for $e = 0.2$ (for $\omega/n \approx 4.5$ and $\omega/n \approx -2.2$), but they are no longer possible for $e = 0.4$.

In Fig. 8b, we show the evolution of the rotation rate for an obliquity $\theta = 45°$, while in Fig. 8c, we show the evolution for $\theta = 90°$. We observe that, as we increase the obliquity, the presence of asynchronous thermal states becomes more difficult. They are still present for $\theta = 45°$ in the circular case ($e = 0$), but for $\theta = 90°$ they are no longer possible for any eccentricity value. On the other hand, for nonzero obliquities, the rotation rate can be trapped in additional spin-orbit resonances, even in the circular case, for $\omega/n = -2, -1, 0, +1, +2$, as it was already reported for Venus (Fig. 2). We hence conclude that capture in the asynchronous thermal states is more likely for low obliquity and small eccentricity, while large obliquities and/or high eccentricities favour the capture in spin-orbit resonances.

### 5.2. Numerical simulations

We now simulate the complete tidal evolution of Kepler-1229 b, by integrating the full set of secular equations of motion for 10 Gyr (Eqs. (15), (16), and (17)). The initial spin state of Earth-like planets is unknown, since large impacts during the formation process can completely modify their rotation period and spin axis orientation (Dones & Tremaine 1993; Kokubo & Ida 2007). In all simulations we assumed that the initial rotation period of Kepler-1229 b is one day, which gives an initial rotation rate of $\omega/n \approx 87$. Since the rotation evolves into an equilibrium state, the choice of its initial value is not critical. For the initial obliquity, we adopted different values between 0° and 180°.

Despite the fact that we take into account the orbital evolution in our model (Eqs. (15) and (17)), we observed that the semi-major axis and the eccentricity remain approximately constant throughout the evolution, which is due to $||G|| \gg ||L||$, where **G** and **L** are the orbital and rotational angular momentum, respectively (Eqs. (1) and (2)). Therefore, we assume for the initial semi-major axis the present value, $a = 0.3125$ au. The eccentricity of Kepler-1229 b is unknown, but this parameter is critical for the spin evolution because it can lead to different final states (Fig. 8). To explore the multiple evolution possibilities, we then adopt three different values for the initial eccentricity, $e = 0, e = 0.2$ and $e = 0.4$.

The results of the numerical simulations for the tidal evolution of Kepler-1229 b over time are shown in Figs. 9 to 11, one for each initial eccentricity value. Since the orbital evolution is negligible for this system, we focus our analysis in the spin evolution of the planet (Eq. (16)). We thus only show the $\omega/n$ ratio (top) and the obliquity (bottom) as a function of time for several different initial obliquities.

<footnote n="10">¹⁰ For $e \neq 0$, we still use Eq. (37), but we get a slightly different mean flux, and thus also slightly different values for $q_0$ and $\sigma_0$.</footnote>
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Secular evolution over time for Kepler 1229 b with $e = 0.0$ and for initial obliquity values between 0° and 180°. We show the $\omega/n$ ratio (top) and the obliquity (bottom). In all simulations, the initial rotation period is 24 h ($\omega/n \approx 87$).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Secular evolution over time for Kepler 1229 b with $e = 0.2$ and for initial obliquity values between 0° and 180°. We show the $\omega/n$ ratio (top) and the obliquity (bottom). In all simulations, the initial rotation period is 24 h ($\omega/n \approx 87$).</figcaption></figure>

In Fig. 9, we show the evolution for $e = 0$ (circular orbit). We observe that, regardless of the initial obliquity, all trajectories end in the prograde asynchronous thermal state in less than about 5 Gyr. More precisely, the spin ends in the prograde state $\omega/n = \omega_s/n + 1 \approx 5.2$ with low obliquity, but not exactly zero ($\theta \approx 14^\circ$). The nonzero final obliquity is somehow surprising, because in the circular case, gravitational and thermal tides usually push the obliquity into either 0° or 180°, respectively (Correia et al. 2003). Of course, in this case the two effects are simultaneously in action, so both torques can cancel for some obliquity (as it happens for the rotation rate).

In Fig. 10, we show the evolution for $e = 0.2$ (small eccentricity). We observe that all simulations end in a stable equilibrium in less than about 4 Gyr. However, contrarily to the circular case, depending on the initial obliquity, in this case the rotation can be either captured in an asynchronous thermal state, or in some spin-orbit resonance. Indeed, in Fig. 8 b, we observe that for $e = 0.2$ and moderate obliquity, the asynchronous thermal states are no longer possible. For initial obliquities close to 180°, the spin evolves into the retrograde asynchronous thermal state $\omega/n \approx 2.2$ with $\theta = 180°$. On the other hand, for initial obliquities lower than 160°, the rotation rate stabilises in the $\omega/n = 4/1$ spin-orbit resonance, while the obliquity settles at $\theta \approx 50°$. Indeed, in Fig. 8 b, we observe that for $e = 0.2$ and an obliquity of 45°, the asynchronous thermal states are no longer possible. We also observe that the first spin-orbit resonance where capture can occur is for $\omega/n = 5/1$. That is why in our simulations the rotation is initially captured in this resonance (Fig. 10). However, in the $5/1$ resonance the obliquity decreases. When the obliquity descends below the threshold $\theta \approx 17°$, the resonant

equilibrium is no longer possible, and the rotation rate decreases. Then, the obliquity increases again and the rotation evolves into the following resonance, $\omega/n = 4/1$, where it can be trapped permanently because the obliquity remains high (Fig. 8 b). As in the circular case, it is remarkable that such high obliquities can be maintained. Valente & Correia (2022) have shown that, when considering gravitational tides alone, some half-integer spin-orbit resonances can excite the obliquity to high values. Here we find that when we additionally include thermal tides, integer resonances can also have the same kind of effect.

In Fig. 11, we show the evolution for $e = 0.4$ (moderate eccentricity). We observe that all simulations end in a stable equilibrium in less than about 1 Gyr. In this case, asynchronous thermal equilibria are not possible for any obliquity value (Fig. 8), so the rotation must evolve into a spin-orbit resonance. Indeed, regardless of the initial obliquity, we observe a succession of resonance captures over time in Fig. 11. For initial obliquities close to 180°, the rotation is directly captured in the $\omega/n = 2/1$ spin-orbit resonance, and the final obliquity evolves into $\theta = 0°$. On the other hand, for initial obliquities lower than 160°, the obliquity first evolves into a high value around $\theta = 75°$, which facilitates the capture in high order spin-orbit resonances, such as $\omega/n = 8/1$. However, as soon as capture occurs, the obliquity decreases and the spin-orbit resonances are exited after the obliquity descends below a certain limit. The last stage for the spin evolution occurs when all trajectories are captured in the resonance $\omega/n = 5/1$, where the final obliquity also evolves into $\theta = 0°$. This spin-orbit resonance was previously unstable at low obliquity for $e = 0.2$, but it is the first resonance to become
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Secular evolution over time for Kepler 1229 b with $e = 0.4$ and for initial obliquity values between 0° and 180°. We show the $\omega/n$ ratio (top) and the obliquity (bottom). In all simulations, the initial rotation period is 24 h ($\omega/n \approx 87$).</figcaption></figure>

stable when $e = 0.4$ (Fig. 8a). We conclude that the spin evolution for $e = 0.4$ is quite similar to the spin evolution in absence of thermal tides (see Valente & Correia 2022). This result is not totally unexpected, because for higher eccentricities the critical semi-major axis is larger (Fig. 5), and for $a < a_c$, thermal tides are dominated by gravitational tides (Fig. 7).

## 6. Discussion and conclusions

In this paper, we have revisited the long-term evolution of Earth-like planets while taking into account both gravitational and thermal atmospheric tides. Preliminary studies have shown that the combined effect of these two kinds of tides may lead to asynchronous equilibrium states, which can help sustain temperate environments and thus more favourable conditions for life (Correia et al. 2008; Cunha et al. 2015; Leconte et al. 2015; Revol et al. 2023). Here, we adopted a vectorial formalism based on Hansen coefficients expansions (Correia & Valente 2022; Valente & Correia 2023), which allowed us to easily follow the spin and the orbital evolution of these planets for any obliquity and orbital eccentricity. This formalism is very general and can be used with a large variety of rheologies and atmospheric models.

Venus is believed to be a unique example in the Solar System of a planet with a balanced equilibrium between gravitational and thermal tides (e.g. Gold & Soter 1969; Dobrovolskis 1980; Correia & Laskar 2001; Auclair-Desrotour et al. 2017a). If this is so, it can therefore be used to gauge some of the unknown parameters in our models. When adopting an Andrade rheology for gravitational tides (Andrade 1910), we found that the time

exponent $\alpha = 0.3$ and a Maxwell relaxation time of $\tau \approx 1500$ yr provide the best fit for the observations. Although these parameters were obtained from dynamical constraints on the present rotation of Venus, they are consistent with those commonly derived for the Earth and rocky planets using alternative methods (e.g. Castillo-Rogez et al. 2011; Tobie et al. 2019; Bolmont et al. 2020).

When we applied our model to Earth-like planets, we confirmed that asynchronous rotation can indeed be expected for planets in the habitable zone of low-mass stars ($0.4 < M_*/M_\odot < 0.85$). Interestingly, we found that Earth-like planets in the habitable zone of stars with masses ~0.8 $M_\odot$ may end up with an equilibrium rotation period of 24 h. However, these planets require semi-major axes around 0.6 au, so the evolution into this equilibrium state is slow; it can take the entire age of the system to achieve depending on the initial rotation rate of the planet. Conversely, if the initial rotation period is already close to 24 h, it will not evolve much. This result is very compelling but needs to be taken with caution because the equilibrium rotations are very sensitive to the models and parameters adopted for the atmosphere. For instance, we did not consider the effect of high-frequency modes, such as the Lamb waves resonance or the effect of varying the thermal forcing profile (e.g. Farhat et al. 2024; Laskar et al. 2024). However, the 24-h equilibrium underlines the importance of studying the combined effect of gravitational and thermal tides in more detail.

We have also shown that the tidal equilibria strongly depend on the orbital eccentricity. In particular, it gives rise to an asymmetry between the sub- and super-synchronous thermal equilibrium states. High eccentricities also enhance the effect of gravitational tides and hence increase the critical semi-major axis that destabilises the asynchronous equilibria provided by thermal tides. On the other hand, the eccentricity increases the number of spin-orbit resonances where the rotation can be trapped, which also helps develop more temperate climates (e.g. Dobrovolskis 2007; Valente & Correia 2022).

Another striking result is that stable nonzero obliquities are a possible outcome of the tidal equilibria. High obliquities promote a more balanced distribution of the stellar flux on the planet surface and thus extend the size of the habitable areas (e.g. Armstrong et al. 2014; Colose et al. 2019; Dobrovolskis 2021). In this paper, we did not explore the obliquity equilibria in great detail because other effects that we did not include in our model, such as core-mantle friction (e.g. Correia et al. 2003; Correia 2006; Cunha et al. 2015) and planetary perturbations (e.g. Laskar & Robutel 1993; Correia & Laskar 2001, 2003; Saillenfest et al. 2019; Su & Lai 2021) can also play an important role in the obliquity evolution. Therefore, these additional effects should be taken into account in future studies on this topic.

The different tidal equilibria for the spin mainly depend on the orbital eccentricity, the rheology of the mantle, and the composition of the atmosphere. The eccentricity is a parameter that can already be accessed for planets detected with the radial-velocity method (e.g. Sagear & Ballard 2023) and with TTVs (e.g. Agol et al. 2021). For the rheology, we expect that Solar System data will provide more constraints on the deformation and dissipation of rocky planets (e.g. Efroimsky & Lainey 2007; Farhat et al. 2022). We also expect that space telescopes such as JWST and ARIEL (Edwards & Tinetti 2022) will be able to provide more insight on the composition of the atmospheres. Inversely, if we are able to measure the spin of a given planet, we can put constraints on the unknown parameters, as we did for Venus (Fig. 1). 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 for the determination of the spin state of these planets.

**Acknowledgements.** We are grateful to the referee Emeline Bolmont for helpful comments. This work was supported by COMPETE 2020 and by FCT – Fundação para a Ciência e a Tecnologia, I.P., Portugal, through the projects SFRH/BD/137958/2018, GRAVITY (PTDC/FIS-AST/7002/2020), PHOBOS (POCI-01-0145-FEDER-029932), ENGAGE SKA (POCI-01-0145-FEDER-022217), and CFisUC (UIDB/04564/2020 and UIDB/04564/2020, with DOI identifiers 10.54499/UIDB/04564/2020 and 10.54499/UIDB/04564/2020, respectively). It was also supported by the French Agence Nationale de la Recherche (AstroMeso ANR-19-CE31-0002-01) and by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Advanced Grant AstroGeo-885250). We acknowledge the Laboratory for Advanced Computing at University of Coimbra (https://www.uc.pt/lca) for providing the resources to perform the numerical simulations.

## References

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

Agol, E., Dorn, C., Grimm, S. L., & et al. 2021, Planet. Sci. J., 2, 1

Anderson, D. L., & Minster, J. B. 1979, Geophys. J. Int., 58, 431

Andrade, E. N. d. C. 1910, Proc. R. Soc. Lond. A, 84, 1

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

Auclair-Desrotour, P., Laskar, J., & Mathis, S. 2017a, A&A, 603, A107

Auclair-Desrotour, P., Laskar, J., Mathis, S., & Correia, A. C. M. 2017b, A&A, 603, A108

Auclair-Desrotour, P., Leconte, J., & Mergny, C. 2019, A&A, 624, A17

Auclair-Desrotour, P., Deitrick, R., & Heng, K. 2022, A&A, 663, A79

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

Bolmont, E., Selsis, F., Owen, J. E., et al. 2017, MNRAS, 464, 3728

Bolmont, E., Breton, S. N., Tobie, G., et al. 2020, A&A, 644, A165

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

Carpenter, R. L. 1964, AJ, 69, 2

Castillo-Rogez, J. C., Efroimsky, M., & Lainey, V. 2011, J. Geophys. Res. (Planets), 116, E09008

Chapman, S., & Lindzen, R. 1970, Atmospheric tides: Thermal and Gravitational (Dordrecht: Reidel)

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

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

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

Correia, A. C. M., & Laskar, J. 2003, Icarus, 163, 24

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

Correia, A. C. M., Laskar, J., & Néron de Surgy, O. 2003, Icarus, 163, 1

Correia, A. C. M., Udry, S., Mayor, M., et al. 2008, A&A, 479, 271

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

Cunha, D., Correia, A. C. M., & Laskar, J. 2015, Int. J. Astrobiol., 14, 233

Cuntz, M., & Wang, Z. 2018, RNAAS, 2, 19

Demory, B.-O., Gillon, M., Madhusudhan, N., & Queloz, D. 2016, MNRAS, 455, 2018

Dobrovolskis, A. R. 1980, Icarus, 41, 18

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

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

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

Dobrovolskis, A. R., & Ingersoll, A. P. 1980, Icarus, 41, 1

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

Dumoulin, C., Tobie, G., Verhoeven, O., Rosenblatt, P., & Rambaux, N. 2017, J. Geophys. Res. (Planets), 122, 1338

Edwards, B., & Tinetti, G. 2022, VizieR Online Data Catalog: J/AJ/164/15

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

Efroimsky, M. 2012b, ApJ, 746, 150

Efroimsky, M., & Lainey, V. 2007, J. Geophys. Res. (Planets), 112, 12003

Farhat, M., Auclair-Desrotour, P., Boué, G., & Laskar, J. 2022, A&A, 665, A1

Farhat, M., Auclair-Desrotour, P., Boué, G., Deitrick, R., & Laskar, J. 2024, A&A, 684, A49

Forget, F., & Leconte, J. 2014, Philos. Trans. Roy. Soc. Lond. Ser. A, 372, 20130084

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

Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, VizieR Online Data Catalog: J/AJ/154/109

Gevorgyan, Y., Boué, G., Ragazzo, C., Ruiz, L. S., & Correia, A. C. M. 2020, Icarus, 343, 113610

Gold, T., & Soter, S. 1969, Icarus, 11, 356

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

Ingersoll, A. P., & Dobrovolskis, A. R. 1978, Nature, 275, 37

Jeffreys, H. 1976, The Earth. Its Origin, History and Physical Constitution

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

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

Koll, D. D. B., & Abbot, D. S. 2016, ApJ, 825, 99

Komacek, T. D., & Abbot, D. S. 2019, ApJ, 871, 245

Kopparapu, R. K., Ramirez, R., Kasting, J. F., et al. 2013, ApJ, 765, 131

Kopparapu, R. K., Wolf, E. T., & Meadows, V. S. 2020, in Planetary Astrobiology, eds. V. S. Meadows, G. N. Arney, B. E. Schmidt, & D. J. Des Marais, 449

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

Laskar, J. 1994, A&A, 287, L9

Laskar, J. 2008, Icarus, 196, 1

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

Laskar, J., Joutel, F., & Boudin, F. 1993, A&A, 270, 522

Laskar, J., Farhat, M., Lantink, M. L., et al. 2024, Sedimentologika, 2, 1271

Leconte, J., Forget, F., Charnay, B., Wordsworth, R., & Pottier, A. 2013, Nature, 504, 268

Leconte, J., Wu, H., Menou, K., & Murray, N. 2015, Science, 347, 632

Louchet, F., & Duval, P. 2009, Int. J. Mater. Res., 100, 1433

Luger, R., & Barnes, R. 2015, Astrobiology, 15, 119

Makarov, V. V. 2012, ApJ, 752, 73

Makarov, V. V., & Efroimsky, M. 2013, ApJ, 764, 27

Margot, J.-L., Campbell, D. B., Giorgini, J. D., et al. 2021, Nat. Astron., 5, 676

Melini, D., Saliby, C., & Spada, G. 2022, Geophys. J. Int., 231, 1502

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

Mocquet, A., Rosenblatt, P., Dehant, V., & Verhoeven, O. 2011, Planet. Space Sci., 59, 1048

Morton, T. D., Bryson, S. T., Coughlin, J. L., et al. 2016, ApJ, 822, 86

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

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

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

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

Revol, A., Bolmont, E., Tobie, G., et al. 2023, A&A, 674, A227

Sabadini, R., Vermeersen, B., & Cambiotti, G. 2016, Global Dynamics of the Earth: Applications of Viscoelastic Relaxation Theory to Solid-Earth and Planetary Geophysics (Dordrecht: Springer Netherlands)

Sagear, S., & Ballard, S. 2023, PNAS, submitted [arXiv:2305.17157]

Saillenfest, M., Laskar, J., & Boué, G. 2019, A&A, 623, A4

Saliby, C., Fienga, A., Briaud, A., Mémin, A., & Herrera, C. 2023, Planet. Space Sci., 231, 105677

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, 195

Siebert, M. 1961, Adv. Geophys., 7, 105

Smith, W. B. 1963, AJ, 68, 15

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

Tobie, G., Grasset, O., Dumoulin, C., & Mocquet, A. 2019, A&A, 630, A70

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

Turcotte, D. L., & Schubert, G. 2002, Geodynamics

Valente, E. F. S., & Correia, A. C. M. 2022, A&A, 665, A130

Valente, E. F. S., & Correia, A. C. M. 2023, A&A, 679, A153

Volland, H. 1990, in Earth's Rotation from Eons to Days, 127

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

Winn, J. N. 2018, in Handbook of Exoplanets, eds. H. J. Deeg, & J. A. Belmonte, 195

Wordsworth, R. 2015, ApJ, 806, 180

Wordsworth, R., & Kreidberg, L. 2022, ARA&A, 60, 159

Wu, H., Murray, N., Menou, K., Lee, C., & Leconte, J. 2023, Sci. Adv., 9, eadd2499

Yang, J., Cowan, N. B., & Abbot, D. S. 2013, ApJ, 771, L45

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

Yoder, C. F. 1997, in Venus II: Geology, Geophysics, Atmosphere, and Solar Wind Environment, 1087