© The Authors 2023

1 Introduction

The five thousand exoplanets discovered so far¹ have revealed a great diversity of worlds. As the number of discoveries continues to grow, an accurate modeling of exoplanets becomes increasingly important. In the context of the arrival of new powerful instruments such as the James Webb Space Telescope (i.e. JWST; Greene et al. 2016) and the Atmospheric Remote-sensing Infrared Exoplanet Large-survey mission (i.e. ARIEL; Tinetti et al. 2021; Edwards & Tinetti 2022) in the characterization of rocky planets, we need to describe the dynamical state of rocky exoplanets with more realistic models by taking their internal structure and their potential atmosphere into account. A large number of the rocky planets discovered so far are in very close-in orbits, and are therefore usually considered to be in a circular and coplanar orbit and with a rotation that is synchronized with their mean motion, showing a permanent dayside. For planets with larger orbits, however, the rotational state and orbital elements (i.e., the semi-major axis, eccentricity, orbital inclination, etc.) evolve on a much longer timescale and are expected to have nonzero eccentricities and/or obliquities. Then, eccentricity or obliquity can trap the spin in a higher rotation state, that is in spin–orbit resonances (hereafter SORs), such as a 3:2 SOR, 2:1 SOR, or higher (e.g., Makarov & Efroimsky 2013; Makarov et al. 2018). If a planet has an atmosphere, another tidal mechanism must be taken into account: the atmospheric thermal tides. These are caused by the differential heating between day- and nightsides (Gold & Soter 1969; Chapman & Lindzen 1970; Dobrovolskis & Ingersoll 1980; Ingersoll & Dobrovolskis 1978; Correia & Laskar 2001; Auclair-Desrotour et al. 2017a). This mechanism is a possible explanation for the current state of Venus as the thermal tides can both desynchronize the planet and increase its obliquity (e.g., Correia & Laskar 2001, 2003a). The rotation rate and the obliquity affect the climate of the planet by influencing the heat distribution. For example, spin rates faster than the synchronization can help prevent atmospheres from collapsing (e.g., Wordsworth 2015) and change the fate of a potential surface water ocean (i.e., complete vaporisation or not; Turbet et al. 2016) through a more effective heat redistribution in the atmosphere. We therefore need a complete dynamical framework with relevant tidal models to determine the rotation states of exoplanets as accurately as possible in the context of future data.

In this article, we use the particular case of Venus to present our recent implementation of planetary tides (Boué & Efroimsky 2019) in a secular code called ESPEM (French acronym for Evolution of Planetary System and Magnetism, Benbakoura et al. 2019; Ahuir et al. 2021). Here we study the case of a Venus-like planet around a Sun-like star. The rotation of Venus is thought to be an equilibrium between the gravitational bodily tides and thermal atmospheric tides (Correia & Laskar 2001, 2003b; Correia et al. 2003). It is therefore a good laboratory for studying these mechanisms. Some unresolved issues still remain, however, such as whether the rotation of Venus is currently in equilibrium, and how it reached its current rotational state. The competition between the tides, gravitational and thermal, strongly depends on the internal state, but in the case of Venus, little is known about its internal structure. This will change with the next incoming mission to Venus, however, as EnVision (Widemann et al. 2020), DAVINCI (Garvin et al. 2022), and VERITAS (Smrekar et al. 2020) will bring valuable data about the internal state of Venus and the thermal atmospheric response of the planet (Bills et al. 2020).

To study the spin evolution of Venus-like planets, and in particular, the capture in SORs, it is necessary to describe the internal structure of the planet and atmosphere well. In particular, Walterová & Běhounková (2020) showed that the internal structure also affects the SOR available by the planet. To ensure a good description in this work, we therefore computed the gravitational bodily tides using the formalism of Kaula (1964) along with the Andrade rheology (Andrade 1910), using rheological parameters constrained from laboratory experiments on olivine (Castillo-Rogez et al. 2011). We computed the thermal tides using the analytical model of Leconte et al. (2015) adapted from the prescription developed by Dobrovolskis & Ingersoll (1980) to reproduce the current state of Venus. We also investigated the effect of the luminosity variation of the star on the equilibrium state between the gravitational and thermal tides.

In Sect. 2, we introduce the tidal model we used for the solid and thermal tides and the implementation in the ESPEM code. In Sect. 3, we discuss the evolution of the spin of a Venus-like planet when we only consider the influence of the solid tide. In particular, we discuss the well-known eccentricity-driven SORs and the less well-known inclination driven SORs. In Sect. 4, we discuss the evolution of the planet taking the thermal tides for the constant and evolving stellar luminosity into account. Finally, we discuss our findings and conclude in Sect. 5.

Fig. 1 Panel a: schematic representation of the atmospheric redistribution caused by the stellar heating on a synchronously rotating planet. The arrows show the movement of the particles of the atmosphere pushed from the hot spot (sub-stellar) toward the cold spots (morning and evening spots). Panel b: delayed deformation of the atmosphere with respect to the position of the star (angle δα in the schema) due to the rotation of the planet. Ω is the spin rate of the planet, and is n the mean motion and the star. Panel c: two tidal contributions, gravitational and thermal, acting in opposition on each other. Figures inspired from Correia & Laskar (2003a).

2 New model of planetary tides in ESPEM

We consider the equilibrium solid tides, which correspond to the mass redistribution of a body (i.e., the planet) under the influence of the gravitational perturbation of a massive (or close-in) orbiting body (i.e. the star). As the planet rotates, the solid bulge will be ahead from the position of the star (as illustrated in red in the Fig. 1c) if the spin of the planet is higher than the mean motion. Then, we consider the so-called thermal tides. These correspond to the mass redistribution of an atmosphere due to stellar heating. In the same manner in which the gradient of the gravitational potential causes the mass redistribution of the body, the thermal tides are raised by the differential heating through the atmosphere (Gold & Soter 1969; Dobrovolskis & Ingersoll 1980; Correia & Laskar 2001). The differential temperature between the day- and nightside causes a pressure gradient and therefore a mass redistribution of the atmosphere. This pressure gradient continuously redistributes the atmospheric particles from the high-temperature side (dayside) to the low-temperature side (nightside). As Fig. 1a shows, the direction of the bulge that forms is parallel to the direction of the heating source (i.e., the star). If the planetary spin is higher than the mean motion, as shown in Fig. 1b, the geometry of the deformation places the atmospheric bulge behind the position of the star by analogy with the solid deformation, while the solid bulge is ahead of the position of the star. The delayed response of the atmosphere caused by its radiative damping affects the dynamics of the planet through viscous coupling at the surface. Figure 1c shows the combination of the gravitational and thermal tides on the planet. The two tides compete until an equilibrium is found.

In the following, we introduce in Sect. 2.1 the formalism we used (Kaula) and in Sect. 2.2 the complex Love number, which allows us to express the tidal potential of the deformed planet in terms of Fourier series. In Sect. 2.3, we extend the notion of the potential Love number to a thermal Love number and detail the way in which we account for the influence of the thermal tides. In Sect. 2.4, we design a homogeneous interior model for the planet that counterbalances the thermal tides for the excitation of Venus. In Sect. 2.5, we discuss the corresponding orbital and rotational equations to finally described the implementation in the ESPEM code (Sect. 2.6).

2.1 Kaula formalism

In order to compute the tidal response of a body to a tidal perturbation, we need to use a formalism that is general enough to encapsulate the frequency-dependent response of a body. This response either requires a decomposition of the tidal potential created by the perturber (hereafter perturbing potential) into Fourier harmonic modes as developed by Kaula (1964), or a time-domain approach as proposed by Correia et al. (2014) and Gevorgyan et al. (2020). Both models allow a study of more complex and realistic rocky and icy bodies (Efroimsky & Makarov 2013; Bolmont et al. 2020a). We used the formalism developed by Darwin (1879) and adapted by Kaula (1964), hereafter the Darwin-Kaula formalism. The Darwin–Kaula theory of bodily tides provides the expression of the perturbing potential of a disturbed body in Fourier series as $$U={\displaystyle \sum _{l=2}^{\infty }{\displaystyle \sum _{m=0}^l{\displaystyle \sum _{p=0}^l{\displaystyle \sum _{q\in \text{ℤ}}{U}_{lmpq}\left(a,e,i,{\sigma }_{lmpq}\right),}}}}$$(1)

with a, e, and i, the semi-major axis, the eccentricity and the inclination respectively. The indexes l, m, p, and q are the indices of the harmonic modes, where l, m are the orders associated with the associated Legendre polynomials, and p, q the order of the Darwin-Kaula Fourier development. Each harmonic mode corresponds to an excitation frequency σ_lmpq, that is, the frequency with which the perturbing potential will affect the deformed body, defined as σ_lmpq = (l − 2 p + q)n − m Ω (with Ω and n the spin rate and the mean motion respectively). Figure 2 shows the contribution of three different modes (l, m, p, q), the (2, 2, 0, 0), (2, 2, 1, 1), and (2, 2, 0, 2) modes, which correspond to the frequencies 2(n − Ω), n − 2Ω, and 4n − 2Ω, respectively. The (l, m, p, q) = (2, 2, 0, 0) mode corresponds to the circular coplanar case (i.e., the semi-diurnal frequency). The (2, 2, 1, 1) and (2, 2, 0, 2) modes are two of the frequencies that are excited when the eccentricity is nonzero.

This formulation is general and fundamental enough to be valid for an arbitrary rheology, and can also be used in the context of the thermal tides (see Sect. 2.3). When the perturber is far enough away, we can keep the development of the gravitational potential at the quadrupolar order only, l = 2 (Makarov & Efroimsky 2013; Mathis & Le Poncin-Lafitte 2009). We restricted the eccentricity expansion numerically up to the order 7, which corresponds to the index 7 in the summation over q and eccentricities up to 0.3.

Fig. 2 Schematic representation of the contribution of the tidal modes l, m, p, q = (2200), (2211), and (2202) in the Kaula formalism. Each bulge represents the tidal deformation under a component of the tidal potential Ulmpq of the perturber (point mass on the right).

2.2 Solid Love number

The response of a planet to the tidal disturbance is quantified using the tidal Love number k₂ (Love 1909). The Love number links the perturbing potential and the additional potential created by the deformed planet in response to the perturbing potential. As we used the quadrupolar component of the tidal potential, we can write each mode as $$U_{2mpq}^{\text{tidal\hspace{0.17em}bulge}}\left({\sigma }_{2mpq}\right)={\overline{k}}_2\left({\sigma }_{2mpq}\right)U_{2mqp}^{\text{tidal\hspace{0.17em}perturbation}}\left({\sigma }_{2mpq}\right),$$(2)

where $U_{2mpq}^{\text{tidal\hspace{0.17em}bulge}}\left({\sigma }_{2mpq}\right)$ corresponds to the quadrupolar component of the potential Eq. (1), (l = 2), and ${\overline{k}}_2\left({\sigma }_{2mpq}\right)$, (hereafter ${\overline{k}}_2$) is related to the amplitude of the complex quadrupolar Love number (Efroimsky 2012a) at the quadrupolar mode l = 2 with a frequency dependence with σ_2mpq.

Physically, the quadrupolar Love number quantifies the response of a body submitted to a periodical external perturbation of frequency σ_2mpq. Here the periodical perturbation corresponds to the tidal potential as described by Eq. (1). It can be written as a complex number, where the real part represents the pure elastic behavior and the imaginary part is the viscous behavior, written as $${\overline{k}}_2=𝔑\left({\overline{k}}_2\right)+i\Im \left({\overline{k}}_2\right)=\left|{\overline{k}}_2\right|\exp \left(-i{ϵ}_2\right).$$(3)

Thus, we can link the phase of the exponent ϵ₂ with the angle between the tidal bulge and the position of the perturber with δ = ϵ₂/2 (Remus et al. 2012). As shown in Fig. 2, each phase is associated with an excitation mode (2mpq).

This complex Love number can be computed for any density, shear modulus and viscosity profiles by integrating the equations of motions and Poisson’s equation relating the displacement, stress, strain and induced potential in the frequency domain assuming a compressible Andrade rheology following the method described in Dumoulin et al. (2017) and Tobie et al. (2019). For a homogeneous solid body, the Love number can be determined from analytical solutions following Efroimsky (2012b).

The tidal torque directly depends on the imaginary part of the Love number. In the circular coplanar case, the tidal torque applied on the planet is expressed as (Kaula 1964; Goldreich 1966; Murray & Dermott 1999) $$T^{\text{grav}}=\frac{3}{2}\frac{𝒢M_{\star }{R}_{\text{p}}^5}{a^6}\Im \left({\overline{k}}_2\right),$$(4)

with 𝒢 the gravitational constant, M_★ the stellar mass, R_p the planetary radius, a the semi-major axis and $\Im \left({\overline{k}}_2\right)$ the imaginary part of the Love number, which can be linked with the well-known dissipation factor Q₂ and the Love number modulus k₂ with $\Im \left({\overline{k}}_2\left(\sigma \right)\right)=-k_2\left(\sigma \right)/Q_2\left(\sigma \right)\text{sign}(\sigma$, Goldreich & Soter 1966; Ogilvie 2014; Bagheri et al. 2022).

Then, we need to model an appropriate Love number ${\overline{k}}_2$ for rocky bodies in order to compute the secular tidal effects. Most models assume that planets are made of weakly viscous fluid (e.g., Hut 1981; Goldreich 1966) even for rocky planets. However, it has been shown that they do not reproduce the correct behavior for highly viscous solid bodies, such as the evolution of their rotation (Henning et al. 2009; Efroimsky & Makarov 2013). We used a more realistic rheological response, the Andrade rheology (Andrade 1910), to better reproduce the behavior of a rocky body under periodical forcing (Castillo-Rogez et al. 2011). The Andrade rheology is an anelastic model built as a combination of dashpots and springs. It is composed with two first components in series, a dashpot and a spring which model the pure viscous damping and the pure elastic rigidity respectively, which correspond to the so-called Maxwell rheology (e.g., Correia et al. 2014). The Maxwell components are linked in series with an infinite number of springs and dashpots in parallel which correspond to the hereditary Andrade property, which retains some aspect of material memory (see Fig. 3 and Efroimsky 2012a for details). This model successfully reproduces a broad range of laboratory measurement of solid behavior under stress and strain, including silicate minerals, metals, and ices (Andrade 1910, 1914; McCarthy & Castillo-Rogez 2013). The rheological profile used in this study was computed with a multilayer model following the method published by Tobie et al. (2005, 2019) and Bolmont et al. (2020b).

The Love number can be computed with the method of Bolmont et al. (2020a). Following Efroimsky (2012b), the complex tidal Love number ${\overline{k}}_2^{\text{grav}}$ is given by $${\overline{k}}_2^{\text{grav}}=\frac{3}{2}\frac{1}{1+A_{2}J/\overline{J}},$$(5)

with J = 1/µ the unrelaxed compliance (with µ the unrelaxed elastic shear modulus in Pa) and A₂ defined by $$A_2=\frac{57J^{-}{}^1}{8\pi 𝒢{\rho }^2{R}_{\text{p}}^2},$$(6)

with ρ the density. $\overline{J}$ is the complex compliance of the material and defined in the formalism of the Andrade rheology with (Castillo-Rogez et al. 2011) $$\overline{J}=J+\beta {\left(i\sigma \right)}^{-\alpha }\text{Γ}\left(1+\alpha \right)-\frac{i}{\eta \sigma },$$(7)

with β a factor that describes the intensity of anelastic friction in the material, Γ the Gamma function, η the shear viscosity, σ the excitation frequency, and α an experimentally fit parameter that which represents the frequency dependence of the transient response. A value of α in the range of 0.23–0.28 allows us to reproduce the dissipation factor and k₂ for the Earth at different frequencies (Tobie et al. 2019).

We studied the case of a Venus-like planet with different end-member temperature profiles. Little is known about the interior structure of Venus. This will likely improve with the upcoming ES EnVision mission (Widemann et al. 2020) and the NASA DaVinci (Garvin et al. 2022) and VERITAS (Smrekar et al. 2020) missions. In the meantime, we considered four possible structures. We used one multilayer profile (referred to as the reference profile) with Earth-like viscosity values as a reference, two other profiles with viscosity values divided by 10 or multiplied by 100 relative to the reference profile and one homogeneous profile (see Sect. 2.4). The multilayer structures can be considered as end members of what we think could be the real interior of Venus (e.g., Bolmont et al. 2020a). The Love numbers associated with the homogeneous profile where computed following the formula described in Bolmont et al. (2020a) and Efroimsky (2012b; see Sect. 2.4). For the multilayer reference structure, we derived the radial density and seismic velocities in the mantle of the planet by using the Perple_X code² (Connolly 2005), which uses a temperature profile from Armann & Tackley (2012) together with the shear modulus profile from the compositional model V1 of Dumoulin et al. (2017). The viscosity was computed as a function of the temperature and pressure profiles as (Dumoulin et al. 2017) $$\eta =\frac{1}{2}{A}_0^{-1}{d}^{2.5}\exp \left(\frac{E_{\text{a}}+P{V}_{\text{a}}}{RT}\right),$$(8)

with E_a and V_a the activation and volume energy, respectively, and A₀ the pre-exponential factor, which are parameters from the Arrhenius equation and depend on the material and d the grain size. The parameters of the dry olivine considered in the upper mantle are E_a = 300 kJ mol⁻¹, A₀ = 6.08 × 10⁻¹⁹ Pa⁻¹ s⁻¹, with a grain size d = 0.68 mm.

Figure 4 shows the internal profiles of the shear modulus µ and the viscosity η for the homogeneous model (see Sect. 2.4), the reference model (hereafter Vref) from Armann & Tackley (2012), and two other models with viscosity profiles obtained by multiplying the viscous reference structure by 0.1 or 100 (denoted V0.1 and V100, respectively).

The metallic core structure was computed using PREM scaled to the Venusian pressure conditions (Dumoulin et al. 2017). The imaginary part of the Love numbers associated with these profiles were computed following Dumoulin et al. (2017) and Bolmont et al. (2020a) and are represented in Fig. 5. A less viscous profile (dash-dotted line Fig. 4) that might correspond to a hotter mantle, will be more dissipative than a more viscous profile (dotted line Fig. 4) for frequencies higher than 10⁻¹¹ s⁻¹ (see Fig. 4).

Fig. 3 Schematic representation of the Andrade anelastic model used in this study (adapted from Renaud & Henning 2018). The two first components in series represent a spring and a dashpot. The elements in parallel represent an infinite number of springs and dashpots.

2.3 Thermal Love number

We considered only the feedback of the tidal bulge of the atmosphere deformed by the pressure gradient at the surface, considering that the atmosphere is perfectly coupled with the surface by viscous friction (e.g., Leconte et al. 2015; Auclair-Desrotour et al. 2019). Thus, we neglected other feedbacks such as the effect of the pressure gradient on the shape of the solid crust and the gravitational anomaly of the atmosphere (see Correia & Laskar 2003a for details). Because the mass redistribution of the atmosphere comes from the surface pressure anomaly, the imaginary part of the complex moment of the surface pressure field $\Im \left(\delta p_{\text{s}}^2\right)$ can be used as a prescription for the imaginary part of the thermal Love number (Leconte et al. 2015; Auclair-Desrotour et al. 2017b). The complex moment of the surface pressure field $\delta p_{\text{s}}^2$ describes the thermal tides amplitude, and therefore, $\Im \left(\delta p_{\text{s}}^2\right)$ can be used to describe the dissipation. Then, we can relate $\Im \left(\delta p_{\text{s}}^2\right)$ to an imaginary thermal Love number $\Im \left(k_2^{\text{thermal}}\right)$. The relation between these two quantities is discussed below.

To calculate the complex moment of the surface pressure field $\delta p_{\text{s}}^2$, we used the work of Leconte et al. (2015), who assumed a Maxwell-like frequency dependence (Ingersoll & Dobrovolskis 1978; Gold & Soter 1969; Auclair-Desrotour et al. 2017b). More realistic frequency dependences have been proposed by Auclair-Desrotour et al. (2019) as a generic formulation and a scaling law, adapted for N₂ atmospheres for different surface pressures. These models will be studied in future developments. Leconte et al. (2015) used a 3D climate model (e.g., Leconte et al. 2013a,b; Forget et al. 2013) specifically tuned for the case of Venus to reproduce the amplitude of the thermal tides on Venus today and used this point to fit an analytical Maxwell-like solution. The analytical formulation of the module of the complex moment of the pressure field is expressed as $$\delta p_{\text{s}}^2=-\frac{q_0}{1+i\frac{\sigma }{2{\omega }_0}},$$(9)

such that the imaginary part can be written as $$\Im \left(\delta p_{\text{s}}^2\right)=\frac{q_0\frac{\sigma }{2{\omega }_0}}{1+{\left(\frac{\sigma }{2{\omega }_0}\right)}^2},$$(10)

with q₀ the amplitude of the quadrupole term of the pressure field at zero frequency, ω₀ the radiative frequency, and σ the excitation frequency. The radiative frequency can be identified with the inverse of the thermal equilibrium timescale. The parameters fit on the GCM simulation of Venus are: q₀ = 201 Pa and ω₀ = 3.77 × 10⁻⁷ s⁻¹ (Leconte et al. 2015). Figure 6 shows the amplitude of the pressure bulge $\left|\delta p_{\text{s}}^2\right|$ as a function of the normalized forcing frequency. The solid curve shows the analytical solution fitted from the values of the amplitude and the phase lag computed with the GCM simulation of Venus of Leconte et al. (2015). As it was not possible to run the specific GCM of Venus for different rotation states, it was not possible to constrain the Maxwell fit better.

The thermal Love number $\Im \left(k_2^{\text{thermal}}\right)$ can be determined from the complex moment of the surface pressure field $\Im \left(\delta p_{\text{s}}^2\right)$ by identification between the thermal and the gravitational torque. The expression of the thermal torque raised by the mass redistribution of the atmosphere is (e.g., Goldreich & Soter 1966; Correia & Laskar 2001, 2003a; Leconte et al. 2015; Auclair-Desrotour et al. 2017a) $$T^{\text{thermal}}=\sqrt{\frac{24\pi }{5}}\frac{M_{\star }}{M_{\text{p}}}\frac{R_{\text{p}}^6}{a^3}\Im \left(\delta p_{\text{s}}^2\right),$$(11)

with M_★ the stellar mass, M_p and R_p the planetary mass and radius respectively, and a the semi-major axis. The thermal equivalent Love number $\Im \left(k_2\right)$ can be written by identification with the expression of the solid torque (Eq. (4)) as $$\Im \left(k_2^{\text{thermal}}\right)=-\sqrt{\frac{32\pi }{15}}\frac{a^3{R}_{\text{p}}}{𝒢M_{\star }{M}_{\text{p}}}\Im \left(\delta p_{\text{s}}^2\right),$$(12)

We note that in contrast to $\Im \left(k_2^{\text{grav}}\right)$, the thermal Love number $\Im \left(k_2^{\text{thermal}}\right)$ depends on the semi-major axis of the planet as the intrinsic response of the atmosphere depends on the flux received by the planet. Then, the two torques, gravitational and thermal, do not have the same dependence as to the semi-major axis, which allows an equilibrium point at which the two tides compensate for each other. The dependence of the mass of the atmosphere is contained in the surface pressure term. We highlight that a more massive atmosphere does not necessarily lead to stronger atmospheric tides. For a more massive atmosphere, the atmospheric layers are more opaque to the stellar flux. Thus, as less stellar flux reaches the surface, the thermal tides are damped for a more massive atmosphere. This effect strongly depends on the composition of the atmosphere and requires a better model to be taken into account. A more massive atmosphere is not investigated in this study. The imaginary Love number $\Im \left(k_2^{\text{thermal}}\right)$ as a function of the excitation frequency for the fitted analytical model (Eq. (12)) is plotted in blue in Fig. 5 in absolute values. In the following, we study the presence equilibrium points as a function of tidal frequency for different internal profiles.

Fig. 4 Shear modulus and viscosity profiles for the multilayer reference model Vreſ and homogeneous model considered here, shown as solid red and black lines, respectively. The dotted and dash-dotted red lines represent the two profiles V0.1 and V100 derived from the multilayer reference model with a viscosity multiplied by ×0.1 and ×100 times, respectively. The viscosity η is computed as in Dumoulin et al. (2017) using Eq. (8) of this work.

Fig. 5 Imaginary part of the gravitational Love number $\Im \left(k_2^{\text{grav}}\right)$ as a function of the excitation frequency of a Venus-like planet for different viscosity profiles. The multilayer reference profile Vreſ derived from Armann & Tackley (2012) is shown as the solid red line. The dotted and dash-dotted red lines represent the V0.1 and V100 profiles derived from the multilayer reference one presented in Fig. 4. In blue we present the imaginary Love number $\Im \left(k_2^{\text{thermal}}\right)$ as a function of the excitation frequency computed with Eq. (12), (in absolute values) associated with the amplitude of the thermal tides presented in Fig. 6. The green curve represents the homogeneous profile described in Sect. 2.4. The vertical dotted black line represents the absolute value of the current frequency state of Venus.

Fig. 6 Amplitude of the pressure bulge $\left|\delta p_{\text{s}}^2\right|$ as a function of the normalized forcing frequency (Ω − n)/Ω. The solid line represents the analytical solution fit for the point of Venus (red dot) computed with Venus GCM simulations (see Leconte et al. 2015). The red bars on the Venus point are not strictly error bars. They represent the dispersion of the pressure bulge at the surface.

2.4 Equilibrium state between gravitational and thermal tides

An equilibrium state between the gravitational and thermal tides can be determined by comparing their imaginary Love numbers. The two tides compensate for each other when the addition of the two imaginary Love numbers is 0 (when the two absolute values are equal; see Fig. 5). Figure 5 shows that the $\Im \left(k_2^{\text{grav}}\right)$ corresponding to the multilayer reference profile (solid red curve) is always higher in amplitude than the $\Im \left(k_2^{\text{thermal}}\right)$, (solid blue curve). Figure 7 shows the spin derivative for the multilayer reference profile, the V0.1 and V100 profiles, and the homogeneous fit profile. Comparing the cases without and with atmosphere (in red and blue, respectively), we find that the solid tides corresponding to the multilayer reference profile are strong enough to compensate for the thermal tides at any frequency. On the one hand, the V0.1 profile is also sufficiently dissipative to compensate for the thermal tides, as it corresponds to a more dissipative interior and thus stronger solid tides. Thus, the spin derivatives associated with the reference profile and the less viscous one are not in equilibrium. The system will therefore evolve to the 1:1 SOR. On other the hand, the V100 profile leads to weaker solid tides. In this case, the spin derivative in Fig. 7 shows that the thermal tides are sufficient to compensate for the gravitational tides, except close to the synchronization, where the solid tides are still strong enough to make the 1:1 SOR stable. This profile, which might correspond to a colder mantle, is not dis-sipative enough to compensate for the thermal tides close to the current state of Venus. This is shown in Fig. 5, where the $\Im \left(k_2^{\text{thermal}}\right)$, (solid blue curve) has a broad range of frequencies (from 10⁻⁷ s⁻¹ to 2 × 10⁻⁵ s⁻¹), where it dominates the gravitational Love number $\Im \left(k_2^{\text{grav}}\right)$ associated with the V100 profile (dashed red curve). Thus, the two intersection points correspond to possible equilibrium states, at which the two tides compensate for each other. The equilibrium points shown in the top panel of Fig. 7 (empty blue circles) correspond to the intersection point at low frequency (about 10⁻⁷ s⁻¹ on Fig. 5). Considering the slope of the derivative, however, this point is not stable. The second point at high frequency (about 4.5 × 10⁻⁵ s⁻¹ on Fig. 5) corresponds to a fast rotation of about three days. This last point is not investigated further because, on the one hand, this state is far from the current state of Venus, and on the other hand, it belongs to the high-frequency regime. In our case, the Maxwelllike frequency-dependent model of thermal tides overestimates the strength of the thermal tides at high frequencies, as the model was fit for the low-frequency regime (Auclair-Desrotour et al. 2019). We therefore consider our approach to be valid in the low-frequency range, that is, for |Ω/n| < 15. As this equilibrium point is very far from the present-day rotation of Venus, it would require a more complex model that is beyond the scope of this study, and we did not investigate it further. A better model, such as the parameterized model proposed by Auclair-Desrotour et al. (2019), will be studied in the future.

None of the profiles reproduce the balance between the two contributions, gravitational and thermal, close to the frequency of Venus Using the method of Bolmont et al. (2020a) described in Sect. 2.1 (Eqs. (5) to (7)), we fit a homogeneous profile (in density, viscosity, and rigidity) that reproduces an equilibrium point at the Venus frequency. As the profile we tried to construct is relatively close to the multilayer profile, we used the parameters given in Table 2 of Bolmont et al. (2020a) for this profile of Venus at α = 0.25 and only fit the value of the viscosity parameter log(η). The other parameters are the rigidity log(µ) = 10.02 in Pa and the ratio of the Andrade and Maxwell time τ_A/τ_M= 0.89 (Castillo-Rogez et al. 2011). The homogeneous profile that fits the thermal Love number $\Im \left(k_2^{\text{thermal}}\right)$ at the Venus frequency is found with log(η) = 22.18 and is plotted in Fig. 5 (green curve). Figure 7 (bottom panel) shows the stable equilibrium point between the gravitational tides associated with the fit profile and the thermal tides (filled blue dot) at the current frequency of Venus as well as two unstable points (empty blue points). These points correspond to the rotation states in which the two tides compensate for each other. We must highlight that the homogeneous profile allows for spin equilibrium at the current frequency of Venus in a very narrow range of internal states. As the interior temperature profile evolves on geologic timescales, it will be relevant to take the associated change of dissipation due to progressive cooling into account, which evolves on timescales of 100 Myr (Bower et al. 2019), or radiogenic decay, tidal heating, and so on, to fully characterize the spin equilibrium. This temperature dependence will be addressed in future studies.

Fig. 7 Spin derivative as a function of the rotation (in terms of Ω/n, Ω and n the planetary spin and mean motion respectively). In the top panel, the red lines correspond to the solid tides, and the blue lines correspond to the cases with solid and atmospheric tides. The solid lines correspond to the reference multilayer profile Vref. The dotted and dash-dotted lines correspond to the V0.1 and V100 profiles, respectively. In the bottom panel, the green line represents the solid tides associated with the homogeneous body (see Sect. 2.2 for details). The blue line corresponds to the cases with the contribution of atmospheric tides. The vertical dashed black line represents the current frequency of Venus in both panels. The dots (filled and empty) represent the equilibrium states between the gravitational and thermal tides (stable and unstable, respectively).

2.5 Secular equations

The secular equations we implemented were derived from the Hamiltonian formalism by Boué & Efroimsky (2019). We used Eqs. (116)–(123) of their work, which were derived within the Darwin-Kaula formalism (see Appendix A). One important hypothesis, that was formulated to derive these equations is the gyroscopic approximation, which implies that the spin rate of a body is much faster than the evolution of the spin-axis orientation. This approximation invalidates the equations when the spin tends to zero for a noncoplanar orbit. A singularity occurs when the spin rate is zero within this approximation (see Boué & Efroimsky 2019). The validity of the equations for an inclined orbit close to the null rotation state will need to be revisited.

In this formalism, the inclination is defined by the angle between the orbital plane and the equatorial plane, which in other words corresponds to the angle between the orbital angular momentum and the planetary spin angular momentum³.

The development of Boué & Efroimsky (2019) also includes the deformation of the secondary under the tidal effect of the primary. We neglected the tidal deformation of the secondary. The resulting secular equations of the spin, eccentricity, and orbital inclination are presented in Appendix A.

2.6 Implementation in the ESPEM code

We implemented the secular equations of Boué & Efroimsky (2019) in the code ESPEM (Benbakoura et al. 2019; Ahuir et al. 2021). This is a secular code integrating the dynamical evolution of a star–planet system. The code takes the coupling between the two layers of low-mass stars into account (convective and radiative layers), as well as the effect of the stellar wind and the torque due to the tides raised by the planet on the convec-tive envelope star and the torque due to the star–planet magnetic interactions for circular and coplanar orbits (Ahuir et al. 2021). The code was only used to compute the angular momentum exchange between the planetary orbit and the stellar radiative core and convective envelope angular momentum (Benbakoura et al. 2019; Ahuir et al. 2021). We have added the tidal torque of the star on the planet within the formalism described in this paper. In addition to the equation for the semi-major axis, we also implemented the equations governing further osculating elements of the planet, such as the eccentricity, orbital inclination, longitude of ascending node, argument of periapsis, and planetary spin. The equations for spin, eccentricity, and inclination, longitude of ascending node and argument of periastron can be found in Eqs. (A.2)–(A.6).

The user needs to provide a data file describing the time evolution of the mass and radius of the star as well as the evolution of the mass and radius of the radiative and convective envelopes, of the moment of inertia and the stellar luminosity. The stellar evolution is computed with evolution files provided by the code STAREVOL (Amard et al. 2016), which gives the internal dissipation and the evolution of the stellar quantities (e.g., the mass, radius of the radiative core, and convective envelope, and luminosity). The evolution of the stellar luminosity is used in Sect. 4.2. The user also needs to specify the initial conditions of the osculating elements of the planet as well as the rotation rate. The code also needs a data file describing the frequency dependence of the real and imaginary parts of the tidal Love number of the planet. This latter file should also provide the mass, radius, and the radius of gyration of the planet (which represents the internal density distribution). The Love numbers provided in these data files were computed with the method described in Sect. 2.4. As explained in Sect. 2.1, several frequencies are excited depending on the eccentricity and inclination. These frequencies were computed for each time step. The real and imaginary parts of the Love number were interpolated linearly from their frequency dependence in the data file. These interpolated values were then used to compute the derivatives of all the quantities mentioned before.

We used the parameters of a Sun-Venus system, that is, a Venus-like mass and radius planet orbiting at 0.723 AU. The parameters we used are listed in Table 1. The initial spin rate is shown from Ω/n = 2.1 for the ESPEM simulations as the SORs we studied are below this spin rate. The longitude of the ascending node and argument of pericenter, were set to zero.

The effect of the stellar tides, stellar wind or the evolution of the stellar layers were not taken into account in this study. We focused here only on the evolution of the rotation state of the planet under the tidal perturbation of a star. We considered this approach appropriate for studying the evolution of the planetary system we consider here. The effect of the Venusian tides inside the Sun can be considered to be negligible as the corresponding evolution timescale is about 10¹⁵ Gyr order of magnitude(e.g., Bolmont & Mathis 2016).

3 Impact of the gravitational tides alone

First, we investigated the secular evolution of the spin, the eccentricity, and the spin inclination of a Venus-like planet orbiting a Sun-like star driven by the gravitational tides alone. In other words, we first neglected the influence of the thermal tides, which is equivalent to first considering an atmosphereless planet. In Sect. 3.1, we discuss spin–orbit resonances for coplanar eccentric orbits, and in Sect. 3.2, we discuss spin–orbit resonances for inclined circular orbits.

3.1 Eccentricity-driven spin-orbit resonances

Hut (1981) showed that if the orbit is eccentric, the planet reaches a pseudo-synchronization state, where the rotation rate of the planet is comparable to the mean motion around the periastron. The use of a model more appropriate for a highly viscous object results in discrete stable spin states in presence of eccentricity, however, in particular, SORs (Makarov & Efroimsky 2013). The higher the eccentricity, the higher the SOR order. A planet beginning its evolution with a high eccentricity and a spin faster than 2.5 times its orbital motion first becomes trapped in the 5:2 SOR. Then, as the eccentricity diminishes, the planet leaves the resonance to be trapped in the lower resonance, the 2:1 SOR, then leaving this configuration for a lower SOR, the 3:2 SOR as the eccentricity continues to decrease. Then, as the eccentricity continues to decrease, the rotation eventually becomes trapped in the 1:1 SOR, also known as the synchronous state, or tidal locking (see also Gomes et al. 2021, with the Creep tidal model).

Figure 8 shows the evolution of the rotation state and the eccentricity of a Venus-like planet with the multilayer internal reference structure (see Sect. 2.2) for three initial eccentricities (0.0, 0.1, and 0.2) in coplanar orbit. The simulations started with an initial semi-major axis of 0.723 AU and an initial rotation period of 100 days. The figure shows that an eccentricity of 0.2 is sufficient to allow the planet to be captured in the 2:1 SOR and the 3:2 SOR for an eccentricity of 0.1. The planet can stay in this SOR as the eccentricity remains high enough throughout the simulation, as shown in the bottom panel of Fig. 8.

The order of the resonance in which the planet capture is shown when we plot the spin derivative as a function of Ω/n (Ω and n the planetary spin and mean motion respectively). Figure 9 shows how this quantity evolves with Ω/n for a fixed eccentricity (Fig. 9a) and how it evolves with Ω/n and for different eccentricities (Fig. 9b). Figure 9a shows that for the circular case (e = 0, dotted blue line), only one value of the spin result in dΩ/dt = 0, and therefore, only one possible equilibrium for the rotation state. The equilibrium is centered at Ω/n = 1, which corresponds to the synchronization state. Higher eccentricities raise other resonances at a higher spin rate. For example, the 0.2 eccentric case (in green on Fig. 9a) shows the 3:2 SOR and the 2:1 SOR in addition to the 1:1 SOR. The 5:2 SOR is also present, but the eccentricity must be higher than 0.25 to keep this configuration stable, as the middle panel shows.

Figure 9b shows the values taken by the spin derivative on a 2D map, as a function of Ω/n and the eccentricity, from 0.0 to 0.3. As we are restricted numerically up to the order 7 in the eccentricity expansion, we computed the evolution up to e = 0.3. This is sufficient because the population of rocky exo-planets does not present extreme eccentricities. The equilibrium points can be found with the null torque in red. The stable equilibrium states must satisfy the condition of a positive torque (in red) to its left and a negative torque (in blue) to its right.

Figure 9b shows that increasing eccentricity allows higher-order SORs. The synchronization is accessible with e = 0, while the 3:2 SOR becomes accessible at e = 0.06, the 2:1 SOR at 0.16, and the 5:2 SOR at 0.26. The eccentric cases of the Fig. 9a are overplotted in the color map with the two horizontal dotted black lines. The evolutions shown in Fig. 8 are plotted in Fig. 9b and 9c with the three colored arrows (with identical colors in the two figures). The spin quickly decreases in the SOR associated with its eccentricity. Because the eccentricity is damped by the tides, the SORs remain until the eccentricity becomes too low to stably maintain these configurations. Figure 9c represents the eccentricity derivative map in the eccentricity versus rotation state plane. In red we show the area in which the eccentricity increases, and in blue the area in which it decreases. In particular, the eccentricity appears to be slightly excited for the e = 0.1 case shown in Fig. 8. This behavior can be explained with the derivative map of Fig. 9c. The eccentricity in the e = 0.2 case of Fig. 8 also appears to be excited before the state when the rotation reached the 2 : 1 SOR.

The timescale of the eccentricity evolution is too long. The departure from resonant states is therefore not shown here.

We confirmed the eccentricity-driven spin-orbit resonances, such as the 1:1, 3:2, 2:1, and 5:2 SORs, and their dependence on the value of the eccentricity. Our results are also consistent with the work of Walterová & Běhounková (2020). We reproduce the eccentricity-driven resonances they showed for the shear modulus and viscosity of our fitted hot profile. Walterová & Běhounková (2020) pointed out that the internal composition of the planet affects the stability of the SORs. Thus, the thermal evolution of the internal structure should be investigated, starting from a warm to a colder profile. As the temperature drives the viscosity and melt fraction of the mantle, the effect of the tidal heating should also be investigated and will be implemented in future developments. Then, the effect of the tidal heating should also be studied, but the tidal dissipation is not thought to be important for the case of Venus. Tidal dissipation is probably stronger for very close-in planets. The effect of the tidal heating of these planets will be the subject of future studies.

Fig. 8 Evolution of a Venus-like planet with an initial rotation of 100 days for three initial eccentricities (i.e., the null eccentricity, and the 0.1 and 0.2 eccentricities). The top panel shows the evolution of the planet rotation rate in terms of Ω/n (Ω and η are the spin and mean motion, respectively). The bottom panel shows the variation in eccentricity Δe (i.e., (e − e(0))/e(0)).

Fig. 9 Spin derivative dΩ/dt(rad s−2). The left panel shows the spin derivative as a function of the rotation state Ω/n (Ω and n are the spin and mean motion, respectively), for different eccentricities. The green dots (filled and empty) represent the equilibrium states, stable (i.e., SORs) and unstable, respectively. The middle and right panels represent the spin derivative and the eccentricity derivative as a function of the rotation state Ω/n and the eccentricity, respectively. The red colored areas depict the positives values, the blue areas depict the negative values, and the red line corresponds to dΩ./dt = 0. The dotted lines represent the two eccentric cases in the left panel (e = 0.1 and e = 0.2). The arrows represent the evolutions presented in Fig. 8.

3.2 Inclination-driven spin–orbit resonances

The inclination-driven SORs has been discussed in Boué et al. (2016) in the context of gas giant planets responding to a Maxwell rheology. We show here that this behavior is also found for rocky planets with a rheology more adapted to rocky planets (Andrade), thus generalizing the findings of Boué et al. (2016). This is the first study of inclination-driven SORs with a realistic rheology for rocky exoplanets that generalizes the first study of Boué et al. (2016) for giant planets, who used a simple Maxwell rheology.

Figure 10 shows the evolution of a Venus-like planet with the multilayer internal reference structure (see Sect. 2.2) after 2 × 10⁸ yr of evolution for three initial inclinations and an initial rotation of 100 days in a circular orbit with the parameters presented in Table 1. For an initial inclination of 5 degrees (blue curve), the tides act to synchronize the rotation of the planet in 5.5 × 10⁷ yr, while the inclination is damped to zero on longer timescales. However, the spin can be trapped in SOR if the initial inclination is high enough. For the initial inclination of 50 and 120 degrees, the planet is captured in the 2:1 SOR for a few 10⁷ yr.

As in the previous section, we investigated how the derivative of the spin depends on rotation state Ω/n and the inclination without eccentricity. Figure 11b shows the spin derivative strength in the plane inclination versus rotation state Ω/n. The equilibrium points can be found with the red curves (null values of the derivative). In the same manner as in Fig. 9, the stable equilibrium states must satisfy the condition of a positive torque (in red) on its left and a negative torque (in blue) on its right for positive values of the rotation state Ω/n, and inversely for negative values of Ω/n. Figure 11a shows that only one equilibrium is possible for low inclinations: synchronous rotation. Increasing the inclination allows other SORs to appear (e.g., the 2:1 SOR). For inclinations higher than 120 degrees, the prograde rotations are no longer equilibrium points, but the retrograde rotations are, such as the −2:1 SOR at Ω/n = −2. A symmetry with respect to a 90-degree inclination exists. This symmetry is clearly visible in the middle panel of Fig. 11b. In particular, the torque at 130-degree inclination is symmetric of the 50-degree inclination. We highlight that no SORs lie above the 2:1 SOR in rotation. Figure 11b shows no SORs close to the Ω/n = 3 or =1.5. Higher spin states were also studied, but as they do not exhibit any SORs. We therefore did not explore a spin rate higher than 100 days. Higher rotation rates require longer timescales to evolve than the present age of the Solar System and are therefore not presented in this paper.

The evolution paths of Fig. 10 are overplotted in Fig. 11 for the four initial inclinations of 5, 50, 120, and 130 degrees. For two of these initial inclinations (50 and 120 degrees), we see a capture in the 2:1 SOR (in Fig. 10 and in Figs. 11b and c). This resonance island is stable for inclinations greater than 15 and lower than 120 degrees. Thus, if the initial spin is higher than Ω/n = 2, the spin is always be damped and trapped in the 2:1 SOR (for an inclination between 15 and 120 degrees). The spin remains in the 2:1 SOR until the inclination becomes too low to stably maintain this configuration.

For an initial inclination of 50 degrees, the inclination appears to be excited by the tides and slightly increases when the spin is higher than the 2:1 SOR. This behavior can be explained with the shape of the inclination derivative di/dt plotted in Fig. 12. It shows a positive-inclination derivative for a spin higher than the 2:1 SOR and an inclination lower than about 80 degrees (bottom right corner of the figure). It also shows that the inclination should also increase when the spin is slightly higher than the synchronization and for an inclination lower than 100 degrees (red area in the vicinity of the 1:1 SOR). This behavior is absent in Fig. 10 because the spin reaches the synchronization very quickly. Higher inclination cases can show an interesting behavior. Figure 11b clearly shows that for a high initial inclination of about 105 degrees, the synchronization state (i.e., Ω/n = 1) is no longer a stable configuration. Then, if the system starts with a positive rotation and a sufficiently high initial inclination, the spin is damped to the antisynchronization state Ω/n = −1, thus a retrograde rotation. As shown in Fig. 12, however, as the spin reaches a negative value, the inclination is driven toward 180 degrees by the tides. This will result in a stable state where the spin is retrograde, in the antisynchronization state, with an inclination of about 180 degrees. This corresponds to a prograde rotation with an orbital inclination of about 0 degrees. This is consistent with the symmetry on the Ω/n = 0 axis and on the 90-degree axis of the derivative map (Fig. 11b).

We investigated spin inclination-driven SORs, such as the 1:1, the 2:1 and their symmetric, the −1:1, and −2:1 SORs, and the evolution of the spin inclination of the planet. We show the range of inclination allowing for the SORs from 0 to 105 degrees for the 1 : 1 and from 20 to 120 degrees for the 2 : 1 SOR. Our simulations in Fig. 10 show the particular behavior of the inclination for a rotation rate above Ω/n = 2, where the inclination appears to be excited by the tides if the inclination is lower than 80 degrees. Finally, the color maps in Fig. 11 show the symmetrical properties of the inclination-driven SORs in Ω/n = 0 and i = 90 degrees.

The effect of the thermal tides of a Venus-like atmosphere for different initial spin inclination is studied in the next section.

Fig. 10 Evolution of a Sun-Venus-like system. The top panel shows the evolution of the rotation state in terms of Ω/n for an initial rotation period of about 100 days. The bottom panel shows the evolution of the orbital inclination for different initial inclinations of about 0, 50, 120, and 130 degrees.

Fig. 11 Spin derivative dΩ/dt (in rad s−2). The left panel shows the spin derivative as a function of the rotation state Ω/n and for different inclinations (5, 50, 120, and 130 degrees). The middle panel represents the value of the spin derivative as a function of Ω/n and the inclination. The red areas depict the positives values, and the blue areas depict the negative values. The dotted black lines of the middle and right panels represent the four inclination values (5, 50, 120, and 130 degrees) plotted in the left panel. The right panel shows a zoom into the square drawn in the middle panel, centered on the 2 : 1 SOR. The dashed black curves depict the paths of the simulations presented in Fig. 10. The gray areas hide the part of the figure close to the null rotation, where the equations used are no longer valid due to the gyroscopic approximation (see Sect. 2.5).

Fig. 12 Same as Fig. 11b, with the inclination derivative di/dt = f(Ω/n, i), in rad s−1, instead of the spin derivative dΩ/dt.

Fig. 13 Panel a: inclination derivative di/dt (in rad/s) as a function of the spin Ω/n (Ω and n are the planetary spin and mean motion respectively) and for different inclinations (from 0 to 180 degrees). Panel b: spin derivative dΩ/dt (in rad/s2) as a function of the spin Ω/n and inclinations (from 0 to 180 degrees). The red lines represent the null derivative in both panels. The arrows show the evolution of the system. They are consistent with the sign of the inclination derivative (panel a) and the spin derivative (panel b). The dotted orange lines represent the null points of the inclination derivative from the panel a (in red in the panel a). The image of Venus at the top of the two panels corresponds to the current state of Venus. The black dot represents the 1:1 synchronization state. The gray area hides the part of the plots close to the null rotation.

4 Venus-like atmospheric tides

As previous studies showed, the current spin state of Venus cannot be reproduced by involving the solid tides alone (Gold & Soter 1969; Dobrovolskis & Ingersoll 1980; Correia & Laskar 2001, 2003a,b; Correia et al. 2003; Leconte et al. 2015). In particular, Correia & Laskar (2001) showed that atmospheric tides can lead to four final rotation states of Venus, one of which is the retrograde rotation observed today. They showed that the current state of Venus cannot be reached for any initial configuration, however. We can consider that the current spin inclination of Venus is either high (about 177.36 degrees) and has a rotation period of 5832.6 h, or a low spin inclination (about 2.64 degrees) and a retrograde rotation. We use the case of Venus as a reference.

The next section (Sect. 4.1) explores the evolution of a Venus-like planet in the spin and inclination parameter space, with a nonevolving atmosphere, a constant luminosity, and a nonevolving internal profile. Section 4.2 explores the effect of the luminosity evolution of the host star, accounting for a simple prescription for the atmospheric evolution.

4.1 Constant luminosity, nonevolving atmosphere

Leconte et al. (2015) fit the parameters of their analytical solution of the pressure bulge (Eq. (9)) to their GCM simulation to model the thermal tides. These parameters are given in Sect. 2.3.

In this section, we investigate the effect of the thermal forcing produced by the host star on the atmosphere. We considered two models for the interior: a multilayer model (introduced in Sect. 2.2) and the fitted homogeneous model (introduced in Sect. 2.4). For the thermal tides, we used the analytical model of thermal tides fit on the present-day Venus (introduced in Sect. 2.3.) The frequency dependence of the corresponding Love numbers is given in Fig. 5.

As discussed in Sect. 2.4, the solid tides corresponding to the multilayer model and its two variants (the V0.1 and V100 profiles) do not allow a stable-equilibrium point close to the current frequency of Venus as they are either too strong or too weak. We fit a homogeneous hot profile, using the method of Bolmont et al. (2020a), in order to find an equilibrium point close to the Venusian frequency (see Sect. 2.4). The shape of the spin derivative in the bottom panel of Fig. 7 shows one stable equilibrium state and two unstable equilibrium states. The two unstable states, close to the synchronization, are also present in the highly viscous V100 profile case. The negative stable spin state was fit to correspond to the retrograde state of Venus. The 1 : 1 synchronous spin state remains stable.

As in Sect. 3.2, we used the derivative maps of dΩ/dt and di/dt as a function of Ω/n and i to represent the evolution of the system in Fig. 13. As discussed in Sect. 2.4, we constrained our study at low spin rates. Because no inclined SORs are higher than the 2:1 (see Sect. 3.2), the initial spin rate of the simulations was set to Ω/n = 2.5 for most of our simulations. We can find the set of initial spins and spin inclinations that can lead to a stable state close to the current rotation state of Venus (with an inclination as high as 180 degrees). Figures 13a and b show the evolution of the system through the map, in which each curve represents an ESPEM simulation. The solid lines show the cases leading to an inclination as high as 180 degrees with a pro-grade rotation, close to the current state of Venus. The dashed lines show the cases leading to the 1 : 1 SOR and an inclination of 0 degrees (i.e., a prograde rotation with a null inclination). The maps show that the atmospheric tides can drive the system toward the high-inclination state and keep the rotation on the pro-grade spin rate (i.e., prograde rotation with a high inclination) if the initial spin inclination is higher than about 150 degrees with a fast initial rotation. This configuration can be reached through the effect of a chaotic motion in the Solar System for the case of Venus (Correia & Laskar 2003a).

Correia & Laskar (2001) argued that the current state of Venus can be described with four final states, depending on the evolution path of the planet. In their work, the paths leading to the current state of Venus either evolved by increasing the spin inclination toward 180 degrees and keeping the spin on a prograde rotation by either decreasing the spin toward retrograde rotation or keeping the spin inclination to zero degrees. In this study, the retrograde rotation can be reached only from the evolution of the spin inclination toward the high-inclination states. None of the paths shown in Fig. 13b crosses the null spin state. Any positive rotation with a low-inclination configuration will drive the system in the synchronous state. As the chaotic effect of a third body will only perturb the spin inclination of the planet it is therefore unlikely that the rotation has crossed the null spin during its evolution given our set of hypotheses.

Fig. 14 Luminosity variation of a Sun-like star over the time from the beginning to the end of the MS. The dash-dotted line represents the simulation start time. The dotted lines represent the current time state.

4.2 Luminosity variation

Dynamical studies of the thermal tides (Correia & Laskar 2001, 2003a; Leconte et al. 2015) have considered a constant luminosity. As the thermal forcing depends on the heat flux of the host star, we also investigated the effect of an evolving luminosity on the rotation evolution of a Venus-like planet. The luminosity evolution of the Sun-like star in ESPEM comes from simulations with the stellar evolution code STAREVOL (Amard et al. 2016). Figure 14 shows the luminosity variation of the Sun-like star we considered. The thermal Love number was computed with the luminosity dependence with the formulation of Auclair-Desrotour et al. (2017b) as $$\Im \left(k_2^{atm}\left(\sigma \right)\right)=-\frac{4}{32}\frac{\kappa \tau \varsigma ϵL_{\star }a}{R_A{T}_0{M}_{\star }R}\frac{\sigma }{{\sigma }^2+{\omega }_0^2},$$(13)

with L_★ and M_★ the stellar mass and luminosity respectively, R the radius of the planet, a the semi-major axis, τ a weight parameter that gives the efficiency of the coupling between the atmosphere and the surface (0 < τ < 1), ς a shape factor depending on the spatial distribution of tidal heat sources, κ the power per mass unit radiated by the atmosphere (where the atmosphere is assumed to behave like a graybody, i.e., Newtonian cooling), ϵ the effective fraction of power absorbed by the atmosphere, σ the excitation frequency, ω₀ the radiative frequency, T₀ the equilibrium surface temperature of the atmosphere, R_A the specific gas constant defined as R_A = R_GP/ℳ_A (R_GP and ℳ_A being the perfect gas constant and the mean molar mass respectively), and α the shape factor depending on the spatial distribution of tidal heat sources. The values of the parameters we used are presented in Table 2.

We started the simulation at 100 Myr regarding the timescale of the rocky planets to form (Chambers 2004). We considered the atmosphere to be fully formed quickly, over the first 1 Myr after the formation of the planet. Then, as the evolution of the atmosphere is very uncertain, we considered it as non-evolving.The following part presents the evolution of the system after 3.6 Gyr of evolution.

Figure 15 shows the spin versus inclination maps with the ESPEM simulation overplotted, in the same manner as the Fig. 13a and 13b. Figures 15a–h show the apparition and the evolution of the equilibrium state close to the current state of Venus (red curve appearing in the top left corner of the map from Figs. 15b–h). In the early stages of the simulations, the stellar flux is lower than today, and gravitational tides dominate the thermal ones. This means that the spin and inclination evolution is mainly driven by the gravitational tides, following a path consistent with Fig. 11c. Figure 15d corresponds to a situation in which an equilibrium close to the current state of Venus is found for the current age of the Solar System. We must emphasize that these maps were set to reproduce the steady state at the current state of Venus for the current solar luminosity. Figures 15e–h show that the Solar luminosity increases faster than the spin state. Thus, as the luminosity increases, the spin never stays in a stable configuration, but continuously evolves toward the stable state. The evolution of the theoretical equilibrium rotation state between the gravitational and thermal tides can be found by finding the rotation rate Ω_eq, which satisfies $k_2^{\text{grav}}\left({\sigma }_{\text{eq}}\right)=k_2^{\text{thermal}}\left({\sigma }_{\text{eq}}\right)$.

Figure 16 shows the evolution of the equilibrium rotation rate in terms of Ω_eq/n and the evolution of the rotation rate Ω/n from the ESPEM simulations shown in Fig. 15. The equilibrium states were determined numerically by finding the rotation state Ω_eq that verifies the equality between the gravitational and the thermal Love number $\Im \left(k_2^{\text{grav}}\left({\text{Ω}}_{eq}/n\right)\right)=\Im \left(k_2^{\text{thermal}}\left({\text{Ω}}_{eq}/n\right)\right)$ over the evolution of the stellar luminosity. We show the evolution of the four simulations that crossed the current spin state of Venus during their evolution in Fig. 15. These cases crossed the equilibrium state close to the current time, but the equilibrium point evolved faster than the rotational state of the simulation. Figure 15h shows that the thermal tides eventually become stronger than the gravitational tides across a large parameter space as the luminosity increases.

In summary, the luminosity evolution leads to two effects. First, the equilibrium changes as the balance between the gravitational tides and thermal tides evolves (thermal tides dominate as the luminosity increases). Second, the planet cannot stay in equilibrium because the timescale of the spin evolution is longer than the timescale of the luminosity evolution. Concerning the first point, as the luminosity increases and the thermal tides becomes stronger, the equilibrium moves to higher spin rates (farther from synchronization). Concerning the second point, the rotation of the planet always chases the equilibria indefinitely (considering a nonevolving atmosphere).

Fig. 15 Same as Fig. 13b for a variable luminosity. Each panel shows the derivation map and ESPEM simulation paths at different time steps. The image of Venus at the top of the plots correspond to the current spin state of Venus. The blue curves correspond to ESPEM simulations evolving with time and luminosity as described in Sect. 4.2.

Fig. 16 Evolution of the rotational equilibrium state Ωeq (in red) evolving with the stellar luminosity (Fig. 14). In blue, we show the four curves corresponding to the rotational evolution from ESPEM that crossed the current spin state of Venus in Fig. 15.

5 Conclusion

We presented the recent implementation of the effect of the tides raised by the star on a telluric Venus-like planet in the code ESPEM. We added the secular evolution of the osculating elements of the planetary orbit (a, e, i, ω, and Ω), that is the semi-major axis, the eccentricity, the inclination, the longitude of ascending node, the argument of periastron, and the planetary spin. We followed the secular equations published by Boué & Efroimsky (2019), which describe the evolution of the osculating elements of the orbit of the planet under tidal perturbations following the Kaula formalism (Kaula 1964). Our implementation includes gravitational and thermal tides, which allowed us to study the tidal effect of an arbitrary atmosphere on an arbitrary planet, provided that the tidal Love numbers k₂ associated with its atmosphere and internal structure are known.

First we focused in Sect. 3.1 on the eccentricity-driven SORs and validate our implementation by finding the 1:1, 3:2, 2:1, and 5:2 SORs, depending on the eccentricity value. Our results are consistent with the findings of Walterová & Bĕhounková (2020). Then we investigated in Sect. 3.2 the inclination-driven SORs, as shown by Boué et al. (2016), here with the Andrade rheology. In particular, we find the 1:1, the 2:1 and their symmetric, the −1:1, and −2:1 SORs.

In Sect. 4.1, we investigated the effect of a thick Venus-like atmosphere with the implementation of the analytical model of Leconte et al. (2015). We used their fit parameters, chosen so that GCM simulations reproduce the current state of Venus. We fit a homogeneous internal structure that allowed gravitational tides to balance thermal tides at the frequency of Venus. We must emphasize that the de-spinning of Venus is a difficult task. Then, we constrained our work to lower initial rotation rates. We find that depending on the initial spin rate and initial spin inclination, either a spin inclination of about zero in the synchronization state or a state close to the current retrograde rotation of Venus results, with a spin rate close to the synchronization and a spin inclination of 180 degrees. The synchronization state (1:1 SOR) is reached when the planet starts with a spin inclination lower than about 120 degrees and a prograde spin in our simulations.

The latter state can be reached when the planet starts either with a high spin inclination (higher than about 120 degrees) and a prograde rotation, or with a spin inclination lower than about 60 degrees and a retrograde rotation. Our results are consistent with the final spin state of Venus found by Correia & Laskar (2001), who computed the evolution of the spin and obliquity of Venus under the solar tides. We point out, however, that Correia & Laskar (2001) used different models for the gravitational tides, which less appropriate than an Andrade rheology for silicate bodies (i.e., the CTL, Hut 1981; Goldreich 1966; Efroimsky & Makarov 2013) and they included the core-mantle friction. The core–mantle friction helps to damp the spin inclination of the planet (Correia & Laskar 2001) and should help the gravitational tides to balance the thermal tides. Furthermore, they also accounted for the chaotic motion in the Solar System (Laskar 1990). In particular, they reported that the chaotic motion helps to transition from low to high inclination. We cannot reproduce this with only two bodies. Further developments of the ESPEM code will include these effects.

In Sect. 4.1, we assumed that the spin state of Venus was in equilibrium state to fit an internal model to the thermal tides. Our results in Sect. 4.2 showed, however, that this may not be the case, and the spin of Venus may still be evolving because of the variation in the solar luminosity. Thus, we investigated the effect of the evolving luminosity on the thermal tides. The evolution of the stellar luminosity leads to a continuous change in the balance between the gravitational and the thermal tides. The rotation of the planet will then continually increase as the luminosity increases. The luminosity evolving faster than the spin, the rotation of the planet will chase the equilibrium state without reaching it. We must highlight that the way in which the thermal tides will continue to increase as the rotation of the planet increase is unclear. In our case, the Maxwell-like frequency-dependent model of thermal tides overestimates the strength of the thermal tides at high frequencies, as the model was fit for the low-frequency regime (Auclair-Desrotour et al. 2019). Further studies with GCM simulations of Venus at higher rotation rates will help to answer this question.

Exploring the complete evolution of the rotation of Venus requires calculating the dynamical evolution of the planet in (at least) three-body simulations. The effect of the evolution of the internal structure and the atmosphere of the planet must also be investigated. The strength of the gravitational and thermal tides, and thus the balance between the two contributions, should have varied strongly during the evolution of the planet since its formation. The current internal structure and atmospheric tides are not sufficiently constrained, however. Future missions to Venus, such as EnVision (Widemann et al. 2020), DAVINCI (Garvin et al. 2022), and VERITAS (Smrekar et al. 2020), will bring valuable data on the internal state of Venus and on the thermal atmospheric response of the planet (Bills et al. 2020). They will help to determine whether the planet is in equilibrium between the gravitational and the thermal tide. These constraints could also help reconstruct the thermal evolution of the planet, which would impact the competition between the gravitational and thermal tides and thus the rotational evolution. Finally, better observations of the atmosphere, together with additional modeling of the Venusian atmosphere, would help constrain the thermal tide. In particular, estimating the response of the atmosphere with a GCM to other frequencies would be extremely helpful.

In the context of exoplanets, we need to consider a relevant model of tides for rocky exoplanets to characterize their surface and potential habitability. We have shown that a relevant tidal model for rocky planet allows a higher spin state than the synchronization, such as eccentricity-driven SORs and also inclination-driven SORs. Planets on a large orbit can keep nonzero eccentricity or obliquity because they evolve on a longer timescale, and can still be trapped in this eccentricity or obliquity-driven SORs. When a planet has an atmosphere, thermal tides can excite the spin inclination to high values because thermal tides drive the spin of Venus in its current state through the chaotic motion of the Solar System. The strength of the thermal tides also depends on the surface pressure, and thus on the total mass of the atmosphere, on the composition that determines the atmospheric absorption, and on the dynamics of the atmosphere. These dependences should be investigated in future studies. We showed that the variation in host star luminosity can also prevent the rotation of a planet from reaching equilibrium between gravitational and thermal tides. This behavior must be further studied for different types of star, that is, different radiation spectra, and with more elaborate models of thermal tides that take the wavelength dependence of the irradiation and the composition of the atmosphere into account. The new generation of instruments, that is, the JWST and ARIEL (Greene et al. 2016; Tinetti et al. 2021; Edwards & Tinetti 2022), will provide valuable data on the atmosphere of rocky worlds. The correct modeling of the dynamical state of exoplanets is then crucial to constrain their surface condition.

Acknowledgements

This work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation under grants 51NF40_182901 and 51NF40_205606. The authors acknowledge the financial support of the SNSF (grant number: 200021_197176). All the members from CEA acknowledge support from GOLF and PLATO CNES grants of the Astrophysics Division at CEA. The computations were performed at University of Geneva on the Baobab and Yggdrasil clusters. This research has made use of NASA’s Astrophysics Data System. The authors thank Drs. Jérémy Leconte, Pierre Auclair-Desrotour and Gwenaël Boué for interesting discussions about the thermal tides.

Appendix A Secular equations

This section presents the secular equations we implemented in the code ESPEM. The secular equations used were developed by Boué & Efroimsky (2019), who revisited the secular equations of Kaula (1964). We used their equations 116 to 123, which describe the secular equations of the osculating elements, a, e, i, $\overline{\text{Ω}}$, ω, and ε the semi-major axis, the eccentricity, the inclination, the longitude of ascending node, the argument of pericenter and the inclination of the spin axis respectively, and θ the spin rate. Here, the inclination i is defined as the angle between the orbital plane and the planet equator. The inclination of the spin axis ε is defined as the inclination of the spin with respect to the inertial frame.

These equations were computed within the gyroscopic approximation, which implies that the spin rate of a body is much faster than the evolution of the spin-axis orientation. This approximation means that considering the limit within which the spin tends to zero cannot be included.

Hereafter, the star is taken as the secondary (subscript_★), and the planet as the primary (subscript _p), C is the inertia momentum, that is, $\text{C=}{M}_{\text{p}}{R}_{\text{p}}^2{r}_g$, with r_g the gyration radius (which represents the internal density distribution). β is the reduced mass β = M_pM_★/(M_p + M_★) and 𝒢 the universal gravitational constant. F_lmp(i) and G_lpq(e) are the inclination and eccentricity polynomials respectively (see Appendix B). The tidal frequency is defined as ${\sigma }_{lmpq}=\left(l-2p+q\right)n-m\dot{\theta }$(with $\dot{\theta }$ and n the spin rate and the mean motion respectively). The phases v_lmpq are defined as $v_{lmpq}=\left(l-2p\right)\omega +\left(l-2p+q\right)ℳ+m\overline{\text{Ω}}$, with ℳ the mean anomaly. k_l is the modulus of the complex Love number $\left|{\overline{k}}_l\right|$ of degree l. The tidal potential energy V₁ is defined from the Hamiltonian formalism as (Boué & Efroimsky 2019) ℋ = ℋ + V₁. The tidal perturbing potential ℛ is computed with ℛ = -V₁/β. The perturbing potential U (Eq. 1) within the Darwin-Kaula formalism can be related to R with $ℛ=-\frac{M_{\star }}{\beta }U$ (Boué & Efroimsky 2019). Then, the tidal perturbing potential ℛ is expressed in the formalism of Kaula (1964) as $$\begin{array}{l}ℛ\left(r,r\prime \right)={{\displaystyle \sum _{l=0}^{+\infty }\left(\frac{R_{\text{p}}}{a\prime }\right)}}^{l+1}\frac{𝒢M\prime }{a\prime }{\left(\frac{R_{\text{p}}}{a\prime }\right)}^l{\displaystyle \sum _{m=0}^l\frac{\left(l-m\right)!}{\left(l+m\right)!}}\left(2-{\delta }_{0,m}\right)\hfill \\ {\displaystyle \sum _{p=0}^l{F}_{lmp}\left(i\prime \right)}{\displaystyle \sum _{q=-\infty }^{\infty }{G}_{lpq}\left(e\prime \right)}{\displaystyle \sum _{h=0}^l{F}_{lmh}\left(i\right)}{\displaystyle \sum _{j=-\infty }^{\infty }{G}_{lhj}\left(e\right)}\hfill \\ k_l\left({\sigma }_{lmpq}\right)\cos \left[\left(v{\prime }_{lmpq}-m\dot{\theta }\prime \right)-\left(v_{lmhj}-m\dot{\theta }\right)-{\epsilon }_l\left({\sigma }_{lmpq}\right)\right],\hfill \end{array}$$(A.1)

where the subscript ′ (i.e., r′, a′, e′, i′, ${\overline{\text{Ω}}}^{\prime }$, and ω′) corresponds to the coordinate of the perturber, and the parameters without subscript (i.e. r, a, e, i, $\overline{\text{Ω}}$, and ω) correspond to the coordinate where the tides are evaluated. In the following, both the perturber and the body on which the tides retro-act are considered to be the star, and therefore we neglect the subscripts. All the following equations come from the Hamiltonian development and the method of Boué & Efroimsky (2019).

The spin derivative equation is expressed from the Hamiltonian formalism as (Boué & Efroimsky 2019) $${\frac{d^2\theta }{d{t}^2}|}_{l=2}=-\frac{\beta }{C}\frac{\partial ℛ}{\partial \overline{\text{Ω}}}$$(A.2)

We carry out the derivative of the perturbing function ℛ, and the secular equation of the spin derivative is written after some algebra as $$\begin{array}{l}{\frac{d^2\theta }{d{t}^2}|}_{l=2}=-\frac{𝒢M_{\star }^2{R}^5}{a^{6}C}{\displaystyle \sum _{m=0}^{2}m}\frac{\left(2-m\right)!}{\left(2+m\right)!}\left(2-{\delta }_{0,m}\right){\displaystyle \sum _{p=0}^2{F}_{2mp}{\left(i\right)}^2}\hfill \\ {\displaystyle \sum _{q=-7}^7{G}_{2pq}{\left(e\right)}^2\Im \left(k_2\left({\sigma }_{2mpq}\right)\right)}\hfill \end{array}$$(A.3)

The eccentricity derivative equation is computed from the Hamiltonian equation as $${\frac{de}{dt}|}_{l=2}=\frac{1-e^2}{n{a}^{2}e}\frac{\partial ℛ}{\partial ℳ}-\frac{\sqrt{1-e^2}}{n{a}^{2}e}\frac{\partial ℛ}{\partial \omega }$$(A.4)

Then, the equation of the eccentricity derivative is computed with the eccentricity squared e² to avoid singularities when the eccentricity tends to zero. Thus, we find $$\begin{array}{l}{\frac{d{e}^2}{dt}|}_{l=2}=-2\sqrt{\left(1-e^2\right)}\sqrt{G\left(M_{\text{p}}+M_{\star }\right)}\frac{M_{\star }}{M_{\text{p}}}\left(\frac{R_{\text{p}}^5}{a^{13/2}}\right)\hfill \\ {\displaystyle \sum _{m=0}^2\frac{\left(2-m\right)!}{\left(2+m\right)!}}\left(2-{\delta }_{0,m}\right){\displaystyle \sum _{p=0}^2{F}_{2mp}{\left(i\right)}^2}\hfill \\ {\displaystyle \sum _{q=-7}^7{G}_{2pq}{\left(e\right)}^2}\left(\sqrt{1-e^2}\left(2-2p+q\right)-\left(2-2p\right)\right)\Im \left(k_2\left({\sigma }_{2mpq}\right)\right)\hfill \end{array}$$(A.5)

The equation of the inclination derivative is also defined from the Hamiltonian formalism as $$\begin{array}{l}{\frac{di}{dt}|}_{l=2}=\frac{\beta }{C\dot{\theta }\sin i}\left(\frac{\partial ℛ}{\partial \omega }-\cos i\frac{\partial ℛ}{\partial \overline{\text{Ω}}}\right)\hfill \\ -\frac{1}{n{a}^2\sqrt{1-e^2}\sin i}\left(\frac{\partial ℛ}{\partial \overline{\text{Ω}}}-\cos i\frac{\partial ℛ}{\partial \omega }\right)\hfill \end{array}$$(A.6)

Then, we carry on the derivative of the perturbing function ℛ, and the equation can be written as $$\begin{array}{l}{\frac{di}{dt}|}_{l=2}=\frac{1}{\sin i}\frac{M_{★}}{M_{\text{p}}}{\left(\frac{R_{\text{p}}}{a}\right)}^5{\displaystyle \sum _{m=0}^2\frac{\left(2-m\right)!}{\left(2+m\right)!}\left(2-{\delta }_{0m}\right)\times }\hfill \\ {\displaystyle \sum _{p=0}^2\left[\frac{\beta n^2{a}^2}{C\dot{\theta }}\left(m\cos i-\left(2-2p\right)\right)-\frac{n}{\sqrt{1-e^2}}\left(\left(2-2p\right)\cos i-m\right)\right]\times }\hfill \\ F_{2mp}{\left(i\right)}^2{\displaystyle \sum _{q=-7}^7{G}_{2pq}{\left(e\right)}^2\Im \left(k_2\left({\sigma }_{2mpq}\right)\right)}\hfill \end{array}$$(A.7)

The equation of the longitude of the ascending node derivative is given as $$\begin{array}{l}{\frac{d\overline{\text{Ω}}}{dt}|}_{l=2}=\left(\frac{\beta \cos i}{C\dot{\theta }\sin i}-\frac{\beta \cos \epsilon \cos \overline{\text{Ω}}}{C\dot{\theta }\sin \epsilon }+\frac{1}{n{a}^2\sqrt{1-e^2}\sin i}\right)\frac{\partial ℛ}{\partial i}\hfill \\ +\frac{\beta \cos \epsilon \sin \overline{\text{Ω}}\cot i}{C\dot{\theta }\sin \epsilon }\frac{\partial ℛ}{\partial \overline{\text{Ω}}}-\frac{\beta \cos \epsilon }{C\dot{\theta }\sin \epsilon }\frac{\sin \overline{\text{Ω}}}{\sin i}\frac{\partial ℛ}{\partial \omega }\hfill \end{array}$$(A.8)

Then, after some algebra, the equation can be written as $$\begin{array}{l}{\frac{d\overline{\text{Ω}}}{dt}|}_{l=2}=\frac{𝒢M_{★}^2{R}_{\text{p}}^5}{a^6}{\displaystyle \sum _{m=0}^2\frac{\left(2-m\right)!}{\left(2+m\right)!}\left(2-{\delta }_{0,m}\right)\{}\hfill \\ \left(\frac{1}{C\dot{\theta }\tan i}-\frac{\cos \overline{\text{Ω}}}{C\dot{\theta }\tan \epsilon }+\frac{1}{\beta n{a}^2\sqrt{1-e^2}\sin i}\right)\times \hfill \\ \frac{1}{2}{\displaystyle \sum _{p=0}^2\frac{\partial F_{2mp}{\left(i\right)}^2}{\partial i}{\displaystyle \sum _{q=-7}^7{G}_{2pq}{\left(e\right)}^2𝔑\left(k_2\left({\sigma }_{2mpq}\right)\right)}}\hfill \\ -\frac{\sin \overline{\text{Ω}}\cot i}{C\dot{\theta }\tan \epsilon }m{\displaystyle \sum _{p=0}^2{F}_2{}_{mp}{\left(i\right)}^2{\displaystyle \sum _{q=-7}^7{G}_{2pq}{\left(e\right)}^2\Im \left(k_2\left({\sigma }_{2mpq}\right)\right)}}\hfill \\ +\frac{1}{C\dot{\theta }\tan \epsilon }\frac{\sin \overline{\text{Ω}}}{\sin i}{\displaystyle \sum _{p=0}^l\left(2-2p\right)F_{2mp}{\left(i\right)}^2{\displaystyle \sum _{q=-7}^7{G}_{2pq}{\left(e\right)}^2\Im \left(k_2\left({\sigma }_{2mpq}\right)\right)}}\}\hfill \end{array}$$(A.9)

The equation of the argument of the periastron derivative is defined as $${\frac{d\omega }{dt}|}_{l=2}=-\frac{\beta }{C\dot{\theta }\sin i}\frac{\partial ℛ}{\partial i}+\frac{\sqrt{1-e^2}}{n{a}^{2}e}\frac{\partial ℛ}{\partial e}-\frac{\cos i}{n{a}^2\sqrt{1-e^2}\sin i}\frac{\partial ℛ}{\partial i}$$(A.10)

Then, the equation can be written as $$\begin{array}{l}{\frac{d\omega }{dt}|}_{l=2}=\frac{𝒢M_{\star }^2{R}_{\text{p}}^5}{a^6}{\displaystyle \sum _{m=0}^2\frac{\left(2-m\right)!}{\left(2+m\right)!}}\frac{\left(2-{\delta }_{0,m}\right)}{2}[\hfill \\ -\left(\frac{1}{C\dot{\theta }\sin i}+\frac{1}{n{a}^2\sqrt{1-e^2}\tan i}\frac{1}{\beta }\right){\displaystyle \sum _{p=0}^2\frac{\partial F_{2mp}{\left(i\right)}^2}{\partial i}}{\displaystyle \sum _{q=-7}^7{G}_{2pq}\left(e\right)}\Re \left(k_2\left({\sigma }_{2mpq}\right)\right)\hfill \\ +\frac{\sqrt{1-e^2}}{n{a}^{2}e}\frac{1}{\beta }{\displaystyle \sum _{p=0}^2{F}_{2mp}\left(i\right)}{\displaystyle \sum _{q=-7}^7\frac{\partial G_{2pq}{\left(e\right)}^2}{\partial e}}\Re \left(k_2\left({\sigma }_{2mpq}\right)\right)]\hfill \end{array}$$(A.11)

Finally, the equation of the inclination of the spin axis is defined as $$\frac{d\epsilon }{dt}=-\frac{\beta }{C\dot{\theta }}\left(\cos \overline{\text{Ω}}\cot i\frac{\partial ℛ}{\partial \overline{\text{Ω}}}+\sin \overline{\text{Ω}}\frac{\partial ℛ}{\partial i}-\frac{\cos \overline{\text{Ω}}}{\sin i}\frac{\partial ℛ}{\partial \omega }\right)$$(A.12)

Then, after some algebra, we find $$\begin{array}{l}\frac{d\epsilon }{dt}=\frac{𝒢M_{★}^2{R}^5}{a^{6}C\dot{\theta }}{\displaystyle \sum _{m=0}^2\frac{\left(2-m\right)!}{\left(2-m\right)!}\left(2-{\delta }_{0,m}\right)(}\hfill \\ m\cos \overline{\text{Ω}}\cot i{\displaystyle \sum _{p=0}^2{F}_{2mp}{\left(i\right)}^2{\displaystyle \sum _{q=-2}^2{G}_{2pq}{\left(e\right)}^2\Im \left(k_2\left({\sigma }_{2mpq}\right)\right)}}\hfill \\ -\frac{1}{2}\sin \overline{\text{Ω}}{\displaystyle \sum _{p=0}^2\frac{\partial F_{2mp}{\left(i\right)}^2}{\partial i}}{\displaystyle \sum _{q=-\infty }^{\infty }{G}_{2pq}{\left(e\right)}^2𝔑\left(k_2\left({\sigma }_{2mpq}\right)\right)}\hfill \\ -\frac{\cos \overline{\text{Ω}}}{\sin i}{\displaystyle \sum _{p=0}^2\left(2-2p\right)F_{2mp}{\left(i\right)}^2{\displaystyle \sum _{q=-2}^2{G}_{2pq}{\left(e\right)}^2\Im \left(k_2\left({\sigma }_{2mpq}\right)\right)}})\hfill \end{array}$$(A.13)

Appendix B Table of the inclination eccentricity from the Fourier development by Kaula, Cayley etc

The inclination and eccentricity polynomials, F_lmp(i) and G_lpq(e), respectively, are given by Eqs. 20, 23, and 24 of Kaula (1961) and are presented in Tables B.1 and B.2. The eccentricity functions are elliptic expansions that can be computed with the Hansen function $X_{l-2p+q}^{-\left(l-1\right),\left(l-2p\right)}$ (Tisserand 1889). These expansions are discussed by Izsak et al. (1964).

We considered eccentricities up to 0.3, which allowed us to consider the eccentricity expansions up to order 7 (see tables Cayley 1861). As the tidal interactions are computed at the quadupolar order l = 2, the index m is constrained between 0 and 2 and p between −7 and 7.

References

All Tables

All Figures

Fig. 1 Panel a: schematic representation of the atmospheric redistribution caused by the stellar heating on a synchronously rotating planet. The arrows show the movement of the particles of the atmosphere pushed from the hot spot (sub-stellar) toward the cold spots (morning and evening spots). Panel b: delayed deformation of the atmosphere with respect to the position of the star (angle δα in the schema) due to the rotation of the planet. Ω is the spin rate of the planet, and is n the mean motion and the star. Panel c: two tidal contributions, gravitational and thermal, acting in opposition on each other. Figures inspired from Correia & Laskar (2003a).

In the text

	Fig. 2 Schematic representation of the contribution of the tidal modes l, m, p, q = (2200), (2211), and (2202) in the Kaula formalism. Each bulge represents the tidal deformation under a component of the tidal potential Ulmpq of the perturber (point mass on the right).
In the text

	Fig. 3 Schematic representation of the Andrade anelastic model used in this study (adapted from Renaud & Henning 2018). The two first components in series represent a spring and a dashpot. The elements in parallel represent an infinite number of springs and dashpots.
In the text

Fig. 4 Shear modulus and viscosity profiles for the multilayer reference model Vreſ and homogeneous model considered here, shown as solid red and black lines, respectively. The dotted and dash-dotted red lines represent the two profiles V0.1 and V100 derived from the multilayer reference model with a viscosity multiplied by ×0.1 and ×100 times, respectively. The viscosity η is computed as in Dumoulin et al. (2017) using Eq. (8) of this work.

In the text

Fig. 5 Imaginary part of the gravitational Love number $\Im \left(k_2^{\text{grav}}\right)$ as a function of the excitation frequency of a Venus-like planet for different viscosity profiles. The multilayer reference profile Vreſ derived from Armann & Tackley (2012) is shown as the solid red line. The dotted and dash-dotted red lines represent the V0.1 and V100 profiles derived from the multilayer reference one presented in Fig. 4. In blue we present the imaginary Love number $\Im \left(k_2^{\text{thermal}}\right)$ as a function of the excitation frequency computed with Eq. (12), (in absolute values) associated with the amplitude of the thermal tides presented in Fig. 6. The green curve represents the homogeneous profile described in Sect. 2.4. The vertical dotted black line represents the absolute value of the current frequency state of Venus.

In the text

Fig. 6 Amplitude of the pressure bulge $\left|\delta p_{\text{s}}^2\right|$ as a function of the normalized forcing frequency (Ω − n)/Ω. The solid line represents the analytical solution fit for the point of Venus (red dot) computed with Venus GCM simulations (see Leconte et al. 2015). The red bars on the Venus point are not strictly error bars. They represent the dispersion of the pressure bulge at the surface.

In the text

Fig. 7 Spin derivative as a function of the rotation (in terms of Ω/n, Ω and n the planetary spin and mean motion respectively). In the top panel, the red lines correspond to the solid tides, and the blue lines correspond to the cases with solid and atmospheric tides. The solid lines correspond to the reference multilayer profile Vref. The dotted and dash-dotted lines correspond to the V0.1 and V100 profiles, respectively. In the bottom panel, the green line represents the solid tides associated with the homogeneous body (see Sect. 2.2 for details). The blue line corresponds to the cases with the contribution of atmospheric tides. The vertical dashed black line represents the current frequency of Venus in both panels. The dots (filled and empty) represent the equilibrium states between the gravitational and thermal tides (stable and unstable, respectively).

In the text

	Fig. 8 Evolution of a Venus-like planet with an initial rotation of 100 days for three initial eccentricities (i.e., the null eccentricity, and the 0.1 and 0.2 eccentricities). The top panel shows the evolution of the planet rotation rate in terms of Ω/n (Ω and η are the spin and mean motion, respectively). The bottom panel shows the variation in eccentricity Δe (i.e., (e − e(0))/e(0)).
In the text

Fig. 9 Spin derivative dΩ/dt(rad s−2). The left panel shows the spin derivative as a function of the rotation state Ω/n (Ω and n are the spin and mean motion, respectively), for different eccentricities. The green dots (filled and empty) represent the equilibrium states, stable (i.e., SORs) and unstable, respectively. The middle and right panels represent the spin derivative and the eccentricity derivative as a function of the rotation state Ω/n and the eccentricity, respectively. The red colored areas depict the positives values, the blue areas depict the negative values, and the red line corresponds to dΩ./dt = 0. The dotted lines represent the two eccentric cases in the left panel (e = 0.1 and e = 0.2). The arrows represent the evolutions presented in Fig. 8.

In the text

	Fig. 10 Evolution of a Sun-Venus-like system. The top panel shows the evolution of the rotation state in terms of Ω/n for an initial rotation period of about 100 days. The bottom panel shows the evolution of the orbital inclination for different initial inclinations of about 0, 50, 120, and 130 degrees.
In the text

Fig. 11 Spin derivative dΩ/dt (in rad s−2). The left panel shows the spin derivative as a function of the rotation state Ω/n and for different inclinations (5, 50, 120, and 130 degrees). The middle panel represents the value of the spin derivative as a function of Ω/n and the inclination. The red areas depict the positives values, and the blue areas depict the negative values. The dotted black lines of the middle and right panels represent the four inclination values (5, 50, 120, and 130 degrees) plotted in the left panel. The right panel shows a zoom into the square drawn in the middle panel, centered on the 2 : 1 SOR. The dashed black curves depict the paths of the simulations presented in Fig. 10. The gray areas hide the part of the figure close to the null rotation, where the equations used are no longer valid due to the gyroscopic approximation (see Sect. 2.5).

In the text

	Fig. 12 Same as Fig. 11b, with the inclination derivative di/dt = f(Ω/n, i), in rad s−1, instead of the spin derivative dΩ/dt.
In the text

Fig. 13 Panel a: inclination derivative di/dt (in rad/s) as a function of the spin Ω/n (Ω and n are the planetary spin and mean motion respectively) and for different inclinations (from 0 to 180 degrees). Panel b: spin derivative dΩ/dt (in rad/s2) as a function of the spin Ω/n and inclinations (from 0 to 180 degrees). The red lines represent the null derivative in both panels. The arrows show the evolution of the system. They are consistent with the sign of the inclination derivative (panel a) and the spin derivative (panel b). The dotted orange lines represent the null points of the inclination derivative from the panel a (in red in the panel a). The image of Venus at the top of the two panels corresponds to the current state of Venus. The black dot represents the 1:1 synchronization state. The gray area hides the part of the plots close to the null rotation.

In the text

	Fig. 14 Luminosity variation of a Sun-like star over the time from the beginning to the end of the MS. The dash-dotted line represents the simulation start time. The dotted lines represent the current time state.
In the text

	Fig. 15 Same as Fig. 13b for a variable luminosity. Each panel shows the derivation map and ESPEM simulation paths at different time steps. The image of Venus at the top of the plots correspond to the current spin state of Venus. The blue curves correspond to ESPEM simulations evolving with time and luminosity as described in Sect. 4.2.
In the text

	Fig. 16 Evolution of the rotational equilibrium state Ωeq (in red) evolving with the stellar luminosity (Fig. 14). In blue, we show the four curves corresponding to the rotational evolution from ESPEM that crossed the current spin state of Venus in Fig. 15.
In the text