© The Authors 2024

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¹. 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 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-1229b, 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 n{a}^2\sqrt{1-e^2}k,\text{\hspace{1em} and \hspace{1em}}E=-\frac{𝒢M_{★}m}{2a},$$(1)

where ß = M_*m(M_* + m)⁻¹ is the reduced mass, a is the semi-major axis, e is the eccentricity, n = [𝒢(M_* + m)a⁻³]^1/2 is the mean motion, 𝒢 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\omega ,$$(2)

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

where r is the distance measured from the planet centre-of-mass, r = ‖r‖ is the norm, $\widehat{r}=r/r$ is the unit vector, r_* is the position of the star, P₂(x) = (3x² − 1)/2 is a Legendre polynomial and ${\widehat{k}}_2$ is the second Love number for potential.

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

$${\widehat{k}}_2(\sigma )=\frac{k_{\text{f}}}{1+\widehat{\mu }(\sigma )},$$(4)

where k_f is the fluid Love number and $\widehat{\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, $\widehat{\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)

$$\widehat{\mu }(\sigma )=\frac{{\tau }_{\upsilon }}{{\tau }_e}{\left[1-\text{i}{\left(\sigma {\tau }_e\right)}^{-1}+{\left(\text{i}\sigma {\tau }_A\right)}^{-\alpha }\text{Γ}(1+\alpha )\right]}^{-1},$$(5)

where τ_v is the viscous relaxation time, τ_e = η/μ is the Maxwell relaxation time, η is the effective viscosity and μ is the rigidity. The parameter α is the time exponent of the Andrade transient creep strain law, which is an empirical adjustable constant that usually ranges from 0.2 ≤ α ≤ 0.4 for rocky planets (Castillo-Rogez et al. 2011). The quantity τ_A is the timescale associated with the Andrade creep and may be termed as the “Andrade” or the “anelastic” time. In general, we have τ_A ≥ τ_e (Makarov 2012), so for simplicity², we assumed that τ_A = τ_e.

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

$${\widehat{k}}_2(\sigma )=a_g(\sigma )-\text{i}{b}_g(\sigma ).$$(6)

This partition is very useful because the imaginary part, b_g(σ), 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_{\text{f}}\frac{ℬ(\sigma )\sigma {\tau }_{\upsilon }}{\mathcal{A}{(\sigma )}^2+ℬ{(\sigma )}^2}=\left(k_{\text{f}}-k_{\text{e}}\right)\frac{ℬ(\sigma )\sigma \tau }{\mathcal{A}{(\sigma )}^2+ℬ{(\sigma )}^2},$$(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)\text{Γ}(1+\alpha )\right],$$(8)

and

$$ℬ(\sigma )=1+|\sigma \tau {|}^{1-\alpha }{\left(\frac{{\tau }_e}{\tau }\right)}^{1-\alpha }\sin \left(\frac{\alpha \pi }{2}\right)\text{Γ}(1+\alpha ),$$(9)

where τ = τ_e + τ_v is the total relaxation time and k_e = k_f τ_e/τ 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=-{\widehat{p}}_2\frac{3}{5\overline{\rho }}{\left(\frac{R}{r}\right)}^3{\left(\frac{a}{r_{★}}\right)}^2{P}_2\left(\widehat{r}\cdot {\widehat{r}}_{★}\right),$$(10)

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

In order to compute ${\widehat{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₀, which gives for the second harmonic of the surface pressure variations (Auclair-Desrotour et al. 2017a, Eq. (166)):

$${\widehat{p}}_2(\sigma )=-\frac{\kappa {\rho }_0{J}_0}{{\sigma }_0+\text{i}\sigma },$$(11)

with κ = 1 − γ⁻¹, where γ is the adiabatic exponent value, ρ₀ is the mean surface density of the atmosphere, and σ₀ 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 ${\widehat{p}}_2$ in its real and imaginary parts as

$${\widehat{p}}_2(\sigma )=a_a(\sigma )-\text{i}{b}_a(\sigma ),$$(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 substellar 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+{\left(\sigma /{\sigma }_0\right)}^2}\sqrt{\frac{10}{3\pi }}{q}_0,$$(13)

where q₀ 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(σ) 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(r)=-{\nabla }_r\left[V_g(r)+V_a(r)\right].$$(14)

The star itself interacts with this field; with mass M_* and located at r = 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_{★}{r}_{★}\times g\left(r_{★}\right),$$(15)

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

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

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

$$\dot{E}=M_{★}{\dot{r}}_{★}\cdot g\left(r_{★}\right).$$(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 canthenaverage 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.

Following Correia & Valente (2022), we let (p, q, s) be a set of orthogonal vectors, such that

$$p=k\times s\text{\hspace{1em} and \hspace{1em}}q=s\times p=k-xs\text{, }$$(18)

with

$$x=\cos \theta =k\cdot s,$$(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 T\rangle =\left(T_p^g+T_p^a\right)p+\left(T_q^g+T_q^a\right)q+\left(T_s^g+T_s^a\right)s,$$(20)

and

$$\langle \dot{E}\rangle =n\left(T_E^g+T_E^a\right),$$(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 T_P 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),

$$\begin{array}{l}{T}_q^g=\frac{3{\mathcal{K}}_g}{32}{\displaystyle \sum _{k=-\infty }^{+\infty }\{3b_g(-kn)\left(1-x^2\right)\left[{\left(X_k^{-3,2}\right)}^2-{\left(X_k^{-3,-2}\right)}^2\right]}\hfill \\ +2b_g(\omega -kn)[{(1+x)}^2(2-x){\left(X_k^{-3,2}\right)}^2\hfill \\ -4x^3{\left(X_k^{-3,0}\right)}^2-{(1-x)}^2(2+x){\left(X_k^{-3,-2}\right)}^2]\hfill \\ +b_g(2\omega -kn)[-4x\left(1-x^2\right){\left(X_k^{-3,0}\right)}^2\hfill \\ +{(1+x)}^3{\left(X_k^{-3,2}\right)}^2-{(1-x)}^3{\left(X_k^{-3,-2}\right)}^2]\},\hfill \end{array}$$(22)

$$\begin{array}{l}{T}_s^g=\frac{3{\mathcal{K}}_g}{32}{\displaystyle \sum _{k=-\infty }^{+\infty }\{2b_g(\omega -kn)(1-x^2)[4x^2{\left(X_k^{-3,0}\right)}^2}\hfill \\ +{(1+x)}^2{\left(X_k^{-3,2}\right)}^2+{(1-x)}^2{\left(X_k^{-3,-2}\right)}^2]\hfill \\ +b_g(2\omega -kn)[4{\left(1-x^2\right)}^2{\left(X_k^{-3,0}\right)}^2\hfill \\ +{(1+x)}^4{\left(X_k^{-3,2}\right)}^2+{(1-x)}^4{\left(X_k^{-3,-2}\right)}^2]\},\hfill \end{array}$$(23)

$$\begin{array}{l}{T}_E^g=\frac{{\mathcal{K}}_g}{64}{\displaystyle \sum _{k=-\infty }^{+\infty }k}\{b_g(-kn)[4{\left(1-3x^2\right)}^2{\left(X_k^{-3,0}\right)}^2\hfill \\ +9{\left(1-x^2\right)}^2\left({\left(X_k^{-3,-2}\right)}^2+{\left(X_k^{-3,2}\right)}^2\right)]\hfill \\ +12b_g(\omega -kn)\left(1-x^2\right)[4x^2{\left(X_k^{-3,0}\right)}^2\hfill \\ +{(1-x)}^2{\left(X_k^{-3,-2}\right)}^2+{(1+x)}^2{\left(X_k^{-3,2}\right)}^2]\hfill \\ +3b_g(2\omega -kn)[4{\left(1-x^2\right)}^2{\left(X_k^{-3,0}\right)}^2\hfill \\ +{(1-x)}^4{\left(X_k^{-3,-2}\right)}^2+{(1+x)}^4{\left(X_k^{-3,2}\right)}^2]\},\hfill \end{array}$$(24)

where

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

and $X_k^{-\ell ,m}(e)$ are the Hansen coefficients, which depend only on the eccentricity³ (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),

$$\begin{array}{l}{T}_q^a=\frac{3{\mathcal{K}}_a}{32}{\displaystyle \sum _{k=-\infty }^{+\infty }\{3b_a(-kn)\left(1-x^2\right)\left[X_k^{-3,2}{X}_k^{-2,2}-X_k^{-3,-2}{X}_k^{-2,-2}\right]}\hfill \\ +2b_a(\omega -kn)[{(1+x)}^2(2-x)X_k^{-3,2}{X}_k^{-2,2}\hfill \\ -4x^3{X}_k^{-3,0}{X}_k^{-2,0}-{(1-x)}^2(2+x)X_k^{-3,-2}{X}_k^{-3,-2}]\hfill \\ +b_a(2\omega -kn)[-4x\left(1-x^2\right)X_k^{-3,0}{X}_k^{-2,0}\hfill \\ +{(1+x)}^3{X}_k^{-3,2}{X}_k^{-2,2}-{(1-x)}^3{X}_k^{-3,-2}{X}_k^{-2,-2}]\},\hfill \end{array}$$(26)

$$\begin{array}{l}{T}_s^a=\frac{3{\mathcal{K}}_a}{32}{\displaystyle \sum _{k=-\infty }^{+\infty }\{2b_a(\omega -kn)(1-x^2)[4x^2{X}_k^{-3,0}{X}_k^{-2,0}}\hfill \\ +{(1+x)}^2{X}_k^{-3,2}{X}_k^{-2,2}+{(1-x)}^2{X}_k^{-3,-2}{X}_k^{-2,-2}]\hfill \\ +b_a(2\omega -kn)[4{\left(1-x^2\right)}^2{X}_k^{-3,0}{X}_k^{-2,0}\hfill \\ +{(1+x)}^4{X}_k^{-3,2}{X}_k^{-2,2}+{(1-x)}^4{X}_k^{-3,-2}{X}_k^{-2,-2}]\},\hfill \end{array}$$(27)

$$\begin{array}{l}{T}_E^a=\frac{{\mathcal{K}}_a}{64}{\displaystyle \sum _{k=-\infty }^{+\infty }k}\{b_a(-kn)[4{\left(1-3x^2\right)}^2{X}_k^{-3,0}{X}_k^{-2,0}\hfill \\ +9{\left(1-x^2\right)}^2\left(X_k^{-3,-2}{X}_k^{-2,-2}+X_k^{-3,2}{X}_k^{-2,2}\right)]\hfill \\ +12b_a(\omega -kn)\left(1-x^2\right)[4x^2{X}_k^{-3,0}{X}_k^{-2,0}\hfill \\ +{(1-x)}^2{X}_k^{-3,-2}{X}_k^{-2,-2}+{(1+x)}^2{X}_k^{-3,2}{X}_k^{-2,2}\hfill \\ +3b_a(2\omega -kn)[4{\left(1-x^2\right)}^2{X}_k^{-3,0}{X}_k^{-2,0}\hfill \\ +{(1-x)}^4{X}_k^{-3,-2}{X}_k^{-2,-2}+{(1+x)}^4{X}_k^{-3,2}{X}_k^{-2,2}]\},\hfill \end{array}$$(28)

where

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

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},$$(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)}}.$$(31)

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

$$a=-\frac{𝒢M_{★}m}{2E},$$(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}}.$$(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 (ω/n = 0.92462) and an obliquity θ = 177.361° (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 σ₀ = 23.8 yr⁻¹ and q₀ = 2.01 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 μ ≈ 80 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 μ ≈ 176 GPa. The average effective viscosity is much more uncertain, η ~ 10²¹ Pas (Karato & Wu 1993), which gives for the Maxwell relaxation time τ_e ~ 200 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 10⁴ yr ago is still ongoing, suggesting that the Earth’s mantle relaxation time can be τ_e ~ 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}.$$(34)

The current eccentricity of Venus (e = 0.007) and its spin state⁴ (ω_e/n = -0.92462, θ_e = 2.639°), are well determined (Margot et al. 2021), where ω_e and θ_e are the current equilibrium rotation and obliquity, respectively. For simplicity, we can consider that Venus evolves in a circular orbit (e ≈ 0) with zero obliquity (θ ≈ 0), and expression (34) simplifies as (Eqs. (23) and (27))

$$\dot{\omega }=-\frac{3{\mathcal{K}}_g}{2C}{b}_g(2\omega -2n)-\frac{3{\mathcal{K}}_a}{2C}{b}_a(2\omega -2n).$$(35)

Assuming that the present rotation of Venus is in an equilibrium state⁵, we must have $\dot{\omega }=0$ for ω = ω_e, and thus

$${\mathcal{K}}_g{b}_g\left(2{\omega }_e-2n\right)=-{\mathcal{K}}_a{b}_a\left(2{\omega }_e-2n\right).$$(36)

The constant parameters 𝒦_g (Eq. (25)) and 𝒦_a (Eq. (29)), along with the surface pressure variation function, b_a (2ω_e − 2n) (Eq. (13)), are well determined for Venus (Table 2). The major source of uncertainty is related to the Love number function, b_g (2ω_e − 2n), which depends on the mantle rheology. For the Andrade model (Sect. 2.1), this uncertainty is transferred to the determination of the α and τ_e parameters (Eq. (7)).

For materials deforming according to an Andrade rheology with linear work hardening, we have α = 1/3 (Louchet & Duval 2009). However, for rocky planets it is commonly accepted that 0.2 < α < 0.4 (Castillo-Rogez et al. 2011). Indeed, for Earth, Tobie et al. (2019) have shown that α ∈ [0.23,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 α = 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 a and τ_e. In Fig. 1, we show the solution of Eq. (36) for these two parameters. We observe that τ_e can vary by several orders of magnitude for the expected range of α ∈ [0.2,0.4]. For α = 0.2 it is about 60 kyr, while for α = 0.4 it is only about 200 yr. Following Melini et al. (2022), we adopt α = 0.3, which gives for the Maxwell relaxation time, τ_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).

Fig. 1 Maxwell relaxation time τe as a function of the a parameter for Venus (Eq. (36)) obtained for the present rotation state.

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 ω = n (synchronous rotation), while the remaining two correspond to the present state and its symmetric, ω − n = ±ω_s, with ω_s/n = 1.92 (Correia & Laskar 2001). The current observed state, ω_e = n – ω_s, is marked with a dot. The unstable equilibria are located at ω − n = ±ω_c, with ω_c/n = 1.31. These unstable points are also important to understand the behaviour of the rotation rate. For ω/n ∈ [−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, ω = n ± ω_s. These results are in agreement with previous studies that used 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. 2c, we show the evolution for θ = 180°. As expected, this case is similar to the case with θ = 0°, since the pair (ω, 180°) behaves like the pair (−ω, 0°) (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 θ = 0°, Venus is expected to evolve into the state ω = n + ω_s, while for θ = 180°, it is expected to end up in the state ω = n− ω_s (present state).

The previous two obliquity values (θ = 0° and 180°) 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 θ = 90°. 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 (ω/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 θ = 90° 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).

Fig. 2 Variation of dω/dt upon ω/n for Venus with e = 0 (Eq. (34)). We show the evolution for different obliquity values: θ=0° (top), θ = 90° (middle), and θ = 180° (bottom plot). The black dots represent the presently observed rotation state.

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 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‖ ≪ ‖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 ≈ 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 (ω/n, θ). We start the integrations with ω/n = 4.5, and different initial obliquity values in the range [0°, 180°]. Tidal effects initially decrease the rotation rate (Eq. (34)), so the evolutionary paths must be followed for decreasing values of ω/n, until they reach an equilibrium position (marked with a black dot).

For the initial low obliquity (θ = 0°), the obliquity does not change much and the rotation rate reaches the asynchronous rotation with ω/n = 2.92, corresponding to the prograde final state (ω = n + ω_s, θ = 0°). For initial high obliquity (θ ≥ 160°), the obliquity evolves into 180° and then the rotation rate reaches the asynchronous rotation with ω/n = 0.92, which is equivalent to the present retrograde final state (ω = n – ω_s, θ = 0°).

For initial obliquities ranging from 20° to 140°, 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 θ ≈ 3°, 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 (ω = n, θ = 0°). This behaviour can be perfectly understood by looking at Fig. 2 a. For near zero obliquities, rotation rates larger than the critical value ω/n = 2.31 move into the asynchronous final state ω = n + ω_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 Venuslike planets using different values for the α and τ_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 (θ ≥ 160°). 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).

Fig. 3 Secular evolution of the spin of Venus in the plane (ω/n, θ) for different initial obliquity values. This plot was obtained with an initial ω/n = 4.5 and by integrating Eq. (16), where the Hansen coefficients $X_k^{\ell ,m}(e)$ are truncated at |k| ≤ 22. The stable equilibrium points are marked with a black dot.

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 Earthlike 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 semimajor 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_* = 1366 W m⁻², they obtain⁶ for the radiative frequency of the atmosphere, σ₀ = 145 yr⁻¹ , and for the amplitude of the atmospheric quadrupole, q₀ = 11.8 mbar. For the gravitational tides we adopt the Andrade model with the parameters α = 0.3 and τ_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_{★}=\langle \frac{L_{★}}{4\pi r_{★}^2}\rangle =\frac{L_{★}}{4\pi a^2\sqrt{1-e^2}}.$$(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_*/M_⊙ < 0.85 as (Cuntz & Wang 2018),

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

with

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

where L_⊙ and M_⊙ 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, ω/n = 1, and zero obliquity, θ = 0° (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)),

$$\begin{array}{l}\dot{\omega }=\frac{3}{2C}{\displaystyle \sum _{k=-\infty }^{+\infty }[{\mathcal{K}}_g{b}_g(2\omega -kn){\left(X_k^{-3,2}\right)}^2}\hfill \\ +{\mathcal{K}}_a{b}_a(2\omega -kn)X_k^{-3,2}{X}_k^{-2,2}].\hfill \end{array}$$(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 ≈ 0.32 au) bellow 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_*/M_⊙ ≈ 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).

Again, three of these points can be stable, one for ω = n (synchronous rotation) and the other two for ω = n ± ω_s. Here, ω_s = |ω − n| is the synodic rotation frequency (that is, the frequency 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 asynchronous 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 𝒦_a/𝒦_g ∝ a³.

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 several 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 eccentricity, 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 ω/n = k/2 (with k ∈ ℤ) arise from the contribution 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, additional 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 outcome 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 spinorbit 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.

Fig. 4 Variation of dω/dt upon ω/n for an Earth-like planet with θ = 0° (Eq. (40)). We show the evolution for different semi-major axes: a/ac = 0.63 (M*/M⊙ ≈ 0.52; top), a = ac = 0.32 au (M*/M⊙ ≈ 0.68; middle), and a/ac = 1.10 (M*/M⊙ ≈ 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).

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 |ω| ≫ n. Then, for zero obliquity, as the rotation approaches the equilibria region (ω ~ n), it first encounters and locks in one of the asynchronous non-resonant equilibria introduced by thermal 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, ω_s/n = |ω/n − 1|, but also the synodic rotation period, P_s = 2π/ω_s (in days). The solid lines follow the state with ω/n > 1 (super-synchronous), while the dashed lines map the state with ω/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 particular 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_⊙ < 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, ω/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. Interestingly, we find that the equilibrium can be close to P_s ≈ 1 day (the current synodic period of the Earth), for some semi-major axes. More precisely, for a ≈ 0.62 au when e = 0 (circular orbit), for a ≈ 0.64 au when e = 0.2 (small eccentricity), for a ≈ 0.66 au when e = 0.3, and for a ≈ 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_* ≈ 0.82 M_⊙. We hence conclude 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 critical 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).

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⁷ constant for each semi-major axis, going from M_*/M_⊙ ≈ 0.7 to 0.8, respectively. We show again the relative synodic frequency, ω_s/n = |ω/n − 1|, together with the corresponding synodic rotation period, P_s = 2π/ω_s (in days). The solid lines show the state with ω/n > 1 (super-synchronous), while the dashed lines give the state with ω/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 ≈ 0.33, which explains why no asynchronous thermal equilibria is observed for e = 0.4 in Fig. 4c.

Fig. 5 Relative synodic frequency, ωs/n = |ω/n − 1| (top), and respective synodic rotation period, Ps = 2π/ω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 ω/n > 1 (super-synchronous), while the dashed lines give the state with ω/n < 1 (sub-synchronous). These plots are shown on a log-scale in the y-axis.

Fig. 6 Relative synodic frequency, ωs/n = |ω/n − 1| (top), and respective synodic rotation period, Ps = 2π/ω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⊙ ≈ 0.70, 0.73, 0.77, 0.81, respectively. The solid lines give the state with ω/n > 1 (super-synchronous), while the dashed lines give the state with ω/n < 1 (sub-synchronous). These plots are shown on a log-scale in the y-axis.

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 axes larger than some critical value, a_c ≈ 0.35 au, which correspond to stars with M_*/M_⊙ ≳ 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 W m⁻², 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₂ atmosphere and p_s = 1 bar), Leconte et al. (2015) also compute the thermal response for the stellar flux F_* = 450 Wm⁻² using general circulation models. They obtain for the radiative frequency of the atmosphere, σ₀ = 74.5 yr⁻¹, and for the amplitude of the atmospheric quadrupole, q₀ = 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_⊙ = 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 W m⁻² and 1366 Wm⁻², one can extrapolate using a straight line that, for Earth-like planets, we have

$$a_c\approx \left(0.4-0.12M_{★}/M_{\odot }\right)\text{ au. }$$(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⁻², while the bottom line corresponds to F_* = 450 Wm⁻². We observe that for M_*/M_⊙ ≳ 0.5 and a ≳ 0.35 au, it is possible to find a wide region where planets can be simultaneously asynchronous and inside the habitable zone⁸.

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 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₀ and c₀ as a function of the semi-major axis, for a dry terrestrial planet with a homogeneous N₂ atmosphere and p_s = 10 bar around a Sun-like star⁹. 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}_{★},$$(42)

which gives

$${\log }_{10}{q}_0=0.345{\log }_{10}{F}_{★}+2.469$$(43)

and

$${\log }_{10}{\sigma }_0=0.430{\log }_{10}{F}_{★}+0.7343,$$(44)

where a is in au, F_* is given in W m⁻², q₀ in Pa, and σ₀ in yr⁻¹. 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 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.

Fig. 7 Critical semi-major axes for Earth-like planets around low-mass stars. We show the results for an atmosphere with ps = 1 bar (in blue) and ps = 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 W m−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).

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 ≈ 2.74 M_⊕, which is able to accrete more gas. Therefore, we assumed a denser N₂ atmosphere for Kepler-1229 b, with p_s = 10 bar. For a circular orbit¹⁰, we get from Eq. (37) that F_* ≈ 606 Wm⁻², which in turn allowed us to obtain for the amplitude of the atmospheric quadrupole, q₀ = 26.9 mbar (Eq. (43)), and for the radiative frequency of the atmosphere, σ₀ = 85.3 yr⁻¹ (Eq. (44)). In Table 2, we summarise the physical and dynamical parameters that we used.

5.1 Equilibrium final states

We first investigate the possible equilibrium final states for the given parameters of Kepler-1229 b. In Fig. 8 a, we show the evolution of the rotation rate for zero obliquity (Eq. (40)) and different values of the eccentricity (e = 0, e = 0.2, 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 ω = n (synchronous rotation), and another two for ω = n ± ω_s (asynchronous thermal states). More precisely, we find ω_s/n ≈ 4.1 (corresponding to P_s ≈ 21 day), which gives one prograde state at ω/n ≈ 5.1, and one retrograde state at ω/n ≈ −3.1. As we increase the eccentricity, additional equilibria appear at spin–orbit resonances ω/n = k/2 (with k ∈ Z), while the asynchronous thermal states tend to disappear. They still exist for e = 0.2 (for ω/n ≈ 4.5 and ω/n ≈ −2.2), but they are no longer possible for e = 0.4.

In Fig. 8 b, we show the evolution of the rotation rate for an obliquity θ = 45°, while in Fig. 8 c, we show the evolution for θ = 90°. We observe that, as we increase the obliquity, the presence of asynchronous thermal states becomes more difficult. They are still present for θ = 45^° in the circular case (e = 0), but for θ = 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 ω/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.

Fig. 8 Variation of dω/dt upon ω/n for Kepler-1229 b. We show the evolution for different obliquity values: θ = 0° (top), θ = 45° (middle), and θ = 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).

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 Earthlike 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-1229b is one day, which gives an initial rotation rate of ω/n ≈ 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‖ ≫ ‖L‖, where G and L are the orbital and rotational angular momentum, respectively (Eqs. (1) and (2)). Therefore, we assume for the initial semi-majoraxis 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 ω/n ratio (top) and the obliquity (bottom) as a function of time for several different initial obliquities.

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 ω/n = ω_s/n + 1 ≈ 5.2 with low obliquity, but not exactly zero (θ ≈ 14°). 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 ω/n ≈ 2.2 with θ = 180°. On the other hand, for initial obliquities lower than 160°, the rotation rate stabilises in the ω/n = 4/1 spin–orbit resonance, while the obliquity settles at θ ≈ 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 ω/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 θ ≈ 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, ω/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 ω/n = 2/1 spin–orbit resonance, and the final obliquity evolves into θ = 0°. On the other hand, for initial obliquities lower than 160°, the obliquity first evolves into a high value around θ = 75°, which facilitates the capture in high order spin–orbit resonances, such as ω/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 ω/n = 5/1, where the final obliquity also evolves into θ = 0°. This spin-orbit resonance was previously unstable at low obliquity for e = 0.2, but it is the first resonance to become 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).

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 ω/n ratio (top) and the obliquity (bottom). In all simulations, the initial rotation period is 24 h (ω/n ≈ 87).

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 ω/n ratio (top) and the obliquity (bottom). In all simulations, the initial rotation period is 24 h (ω/n ≈ 87).

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 ω/n ratio (top) and the obliquity (bottom). In all simulations, the initial rotation period is 24 h (ω/n ≈ 87).

6 Discussion and conclusions

In this paper, we have revisited the long-term evolution of Earthlike 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 α = 0.3 and a Maxwell relaxation time of τ ≈ 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_⊙ < 0. 85). Interestingly, we found that Earth-like planets in the habitable zone of stars with masses ~0.8 M_⊙ 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 radialvelocity 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 UIDP/04564/2020, with DOI identifiers 10.54499/UIDB/04564/2020 and 10.54499/UIDP/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

All Tables

All Figures

	Fig. 1 Maxwell relaxation time τe as a function of the a parameter for Venus (Eq. (36)) obtained for the present rotation state.
In the text

	Fig. 2 Variation of dω/dt upon ω/n for Venus with e = 0 (Eq. (34)). We show the evolution for different obliquity values: θ=0° (top), θ = 90° (middle), and θ = 180° (bottom plot). The black dots represent the presently observed rotation state.
In the text

	Fig. 3 Secular evolution of the spin of Venus in the plane (ω/n, θ) for different initial obliquity values. This plot was obtained with an initial ω/n = 4.5 and by integrating Eq. (16), where the Hansen coefficients $X_k^{\ell ,m}(e)$ are truncated at |k| ≤ 22. The stable equilibrium points are marked with a black dot.
In the text

	Fig. 4 Variation of dω/dt upon ω/n for an Earth-like planet with θ = 0° (Eq. (40)). We show the evolution for different semi-major axes: a/ac = 0.63 (M*/M⊙ ≈ 0.52; top), a = ac = 0.32 au (M*/M⊙ ≈ 0.68; middle), and a/ac = 1.10 (M*/M⊙ ≈ 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).
In the text

Fig. 5 Relative synodic frequency, ωs/n = |ω/n − 1| (top), and respective synodic rotation period, Ps = 2π/ω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 ω/n > 1 (super-synchronous), while the dashed lines give the state with ω/n < 1 (sub-synchronous). These plots are shown on a log-scale in the y-axis.

In the text

Fig. 6 Relative synodic frequency, ωs/n = |ω/n − 1| (top), and respective synodic rotation period, Ps = 2π/ω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⊙ ≈ 0.70, 0.73, 0.77, 0.81, respectively. The solid lines give the state with ω/n > 1 (super-synchronous), while the dashed lines give the state with ω/n < 1 (sub-synchronous). These plots are shown on a log-scale in the y-axis.

In the text

Fig. 7 Critical semi-major axes for Earth-like planets around low-mass stars. We show the results for an atmosphere with ps = 1 bar (in blue) and ps = 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 W m−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).

In the text

	Fig. 8 Variation of dω/dt upon ω/n for Kepler-1229 b. We show the evolution for different obliquity values: θ = 0° (top), θ = 45° (middle), and θ = 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).
In the text

	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 ω/n ratio (top) and the obliquity (bottom). In all simulations, the initial rotation period is 24 h (ω/n ≈ 87).
In the text

	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 ω/n ratio (top) and the obliquity (bottom). In all simulations, the initial rotation period is 24 h (ω/n ≈ 87).
In the text

	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 ω/n ratio (top) and the obliquity (bottom). In all simulations, the initial rotation period is 24 h (ω/n ≈ 87).
In the text