Issue 
A&A
Volume 674, June 2023



Article Number  A227  
Number of page(s)  18  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202245790  
Published online  27 June 2023 
Spin evolution of Venuslike planets subjected to gravitational and thermal tides
^{1}
Observatoire de Genève, Centre pour la Vie dans l’Univers, Université de Genève,
51 Chemin Pegasis,
1290
Sauverny, Switzerland
email: alexandre.revol@unige.ch; emeline.bolmont@unige.ch
^{2}
Laboratoire de Planétologie et Géosciences, UMRCNRS 6112, Nantes Université,
2 rue de la Houssinière,
BP 92208,
44322
Nantes Cedex 3, France
^{3}
AIM, CEA, CNRS, Université ParisSaclay, Université Paris Diderot,
Sorbonne Paris Cité,
91191
GifsurYvette Cedex, France
Received:
20
December
2022
Accepted:
28
February
2023
Context. The arrival of powerful instruments will provide valuable data for the characterization of rocky exoplanets. Rocky planets are mostly found in closein orbits. They are therefore usually close to the circularcoplanar orbital state and are thus considered to be in a tidally locked synchronous spin state. For planets with larger orbits, however, exoplanets should still have nonzero eccentricities and/or obliquities, and realistic models of tides for rocky planets can allow for higher spin states than the synchronization state in the presence of eccentricities or obliquities.
Aims. This work explores the secular evolution of a star–planet system under tidal interactions, both gravitational and thermal, induced by the quadrupolar component of the gravitational potential and the irradiation of the planetary surface, respectively. We show the possible spin–orbit evolution and resonances for eccentric orbits and explore the possibility of spinorbit resonances raised by the obliquity of the planet. Then, we focus on the additional effect of a thick atmosphere on the possible resulting spin equilibrium states and explore the effect of the evolution of the stellar luminosity.
Methods. We implemented the general secular evolution equations of tidal interactions in the secular code called ESPEM. In particular, we focus here on the tides raised by a star on a rocky planet and consider the effect of the presence of an atmosphere, neglecting the contribution of the stellar tide. The solid part of the tides was modeled with an anelastic rheology (Andrade model), while the atmospheric tides were modeled with an analytical formulation that was fit using a global climate model simulation. We focused on a SunVenuslike system in terms of stellar parameters, orbital configuration and planet size and mass. The SunVenus system is a good laboratory for studying and comparing the possible effect of atmospheric tides, and thus to explore the possible spin state of potential Venuslike exoplanets.
Results. The formalism of Kaula associated with an Andrade rheology allows spin orbit resonances on pure rocky worlds. Similarly to the highorder spin–orbit resonances induced by eccentricity, the spin obliquity allows the excitation of highfrequency Fourier modes that allow some spinorbit resonances to be stable. If the planet has a dense atmosphere, like that of Venus, another mechanism, the thermal tides, can counterbalance the effect of the gravitational tides. We found that thermal tides change the evolution of the spin of the planet, including the capture in spin–orbit resonances. If the spin inclination is high enough, thermal tides can drive the spin toward an antisynchronization state, that is, a the 1:1 spin–orbit resonance with an obliquity of 180 degrees.
Conclusions. Through our improvement of the gravitational and thermal tidal models, we can determine the dynamical state of exoplanets better, especially if they hold a thick atmosphere. In particular, the contribution of the atmospheric tides allows us to reproduce the spin state of Venus at a constant stellar luminosity. Our simulations have shown that the secular evolution of the spin and obliquity can lead to a retrograde spin of the Venuslike planet if the system starts from a high spin obliquity, in agreement with previous studies. The perturbing effect of a third body is still needed to determine the current state of Venus starting from a low initial obliquity. When the luminosity evolution of the Sun is taken into account, the picture changes. We find that the planet never reaches equilibrium: the timescale of the rotation evolution is longer than the luminosity variation timescale, which suggests that Venus may never reach a spin equilibrium state, but may still evolve.
Key words: planets and satellites: terrestrial planets / planets and satellites: dynamical evolution and stability / planet–star interactions
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1 Introduction
The five thousand exoplanets discovered so far^{1} 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 Remotesensing Infrared Exoplanet Largesurvey 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 closein 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 semimajor 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; AuclairDesrotour 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 Venuslike planet around a Sunlike 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 Venuslike 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 (CastilloRogez 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 Venuslike planet when we only consider the influence of the solid tide. In particular, we discuss the wellknown eccentricitydriven SORs and the less wellknown 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 (substellar) 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 closein) 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 socalled 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 hightemperature side (dayside) to the lowtemperature 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 frequencydependent 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 timedomain 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 DarwinKaula formalism. The Darwin–Kaula theory of bodily tides provides the expression of the perturbing potential of a disturbed body in Fourier series as (1)
with a, e, and i, the semimajor 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 DarwinKaula 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 semidiurnal 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 PoncinLafitte 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 U_{lmpq} 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_{2} (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 (2)
where corresponds to the quadrupolar component of the potential Eq. (1), (l = 2), and , (hereafter ) 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 (3)
Thus, we can link the phase of the exponent ϵ_{2} with the angle between the tidal bulge and the position of the perturber with δ = ϵ_{2}/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) (4)
with 𝒢 the gravitational constant, M_{★} the stellar mass, R_{p} the planetary radius, a the semimajor axis and the imaginary part of the Love number, which can be linked with the wellknown dissipation factor Q_{2} and the Love number modulus k_{2} with , Goldreich & Soter 1966; Ogilvie 2014; Bagheri et al. 2022).
Then, we need to model an appropriate Love number 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 (CastilloRogez 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 socalled 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 & CastilloRogez 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 is given by (5)
with J = 1/µ the unrelaxed compliance (with µ the unrelaxed elastic shear modulus in Pa) and A_{2} defined by (6)
with ρ the density. is the complex compliance of the material and defined in the formalism of the Andrade rheology with (CastilloRogez et al. 2011) (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_{2} for the Earth at different frequencies (Tobie et al. 2019).
We studied the case of a Venuslike planet with different endmember 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 Earthlike 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^{2} (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) (8)
with E_{a} and V_{a} the activation and volume energy, respectively, and A_{0} the preexponential 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^{−1}, A_{0} = 6.08 × 10^{−19} Pa^{−1} s^{−1}, 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 (dashdotted 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^{−11} s^{−1} (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; AuclairDesrotour 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 can be used as a prescription for the imaginary part of the thermal Love number (Leconte et al. 2015; AuclairDesrotour et al. 2017b). The complex moment of the surface pressure field describes the thermal tides amplitude, and therefore, can be used to describe the dissipation. Then, we can relate to an imaginary thermal Love number . The relation between these two quantities is discussed below.
To calculate the complex moment of the surface pressure field , we used the work of Leconte et al. (2015), who assumed a Maxwelllike frequency dependence (Ingersoll & Dobrovolskis 1978; Gold & Soter 1969; AuclairDesrotour et al. 2017b). More realistic frequency dependences have been proposed by AuclairDesrotour et al. (2019) as a generic formulation and a scaling law, adapted for N_{2} 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 Maxwelllike solution. The analytical formulation of the module of the complex moment of the pressure field is expressed as (9)
such that the imaginary part can be written as (10)
with q_{0} the amplitude of the quadrupole term of the pressure field at zero frequency, ω_{0} 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_{0} = 201 Pa and ω_{0} = 3.77 × 10^{−7} s^{−1} (Leconte et al. 2015). Figure 6 shows the amplitude of the pressure bulge 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 can be determined from the complex moment of the surface pressure field 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; AuclairDesrotour et al. 2017a) (11)
with M_{★} the stellar mass, M_{p} and R_{p} the planetary mass and radius respectively, and a the semimajor axis. The thermal equivalent Love number can be written by identification with the expression of the solid torque (Eq. (4)) as (12)
We note that in contrast to , the thermal Love number depends on the semimajor 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 semimajor 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 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 dashdotted 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 as a function of the excitation frequency of a Venuslike 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 dashdotted 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 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 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 corresponding to the multilayer reference profile (solid red curve) is always higher in amplitude than the , (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 dissipative enough to compensate for the thermal tides close to the current state of Venus. This is shown in Fig. 5, where the , (solid blue curve) has a broad range of frequencies (from 10^{−7} s^{−1} to 2 × 10^{−5} s^{−1}), where it dominates the gravitational Love number 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^{−7} s^{−1} 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^{−5} s^{−1} 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 highfrequency regime. In our case, the Maxwelllike frequencydependent model of thermal tides overestimates the strength of the thermal tides at high frequencies, as the model was fit for the lowfrequency regime (AuclairDesrotour et al. 2019). We therefore consider our approach to be valid in the lowfrequency range, that is, for Ω/n < 15. As this equilibrium point is very far from the presentday 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 AuclairDesrotour 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 (CastilloRogez et al. 2011). The homogeneous profile that fits the thermal Love number 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 dashdotted 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 DarwinKaula 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 spinaxis 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^{3}.
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 lowmass 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 convective 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 semimajor 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 SunVenus system, that is, a Venuslike 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^{15} Gyr order of magnitude(e.g., Bolmont & Mathis 2016).
Numerical values used in the case of a SunVenuslike system.
3 Impact of the gravitational tides alone
First, we investigated the secular evolution of the spin, the eccentricity, and the spin inclination of a Venuslike planet orbiting a Sunlike 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 Eccentricitydriven spinorbit resonances
Hut (1981) showed that if the orbit is eccentric, the planet reaches a pseudosynchronization 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 Venuslike 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 semimajor 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 exoplanets 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 higherorder 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 eccentricitydriven spinorbit 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 eccentricitydriven 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 closein planets. The effect of the tidal heating of these planets will be the subject of future studies.
Fig. 8 Evolution of a Venuslike 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 Inclinationdriven spin–orbit resonances
The inclinationdriven 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 inclinationdriven 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 Venuslike planet with the multilayer internal reference structure (see Sect. 2.2) after 2 × 10^{8} 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^{7} 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^{7} 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 90degree inclination exists. This symmetry is clearly visible in the middle panel of Fig. 11b. In particular, the torque at 130degree inclination is symmetric of the 50degree 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 positiveinclination 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 90degree axis of the derivative map (Fig. 11b).
We investigated spin inclinationdriven 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 inclinationdriven SORs in Ω/n = 0 and i = 90 degrees.
The effect of the thermal tides of a Venuslike atmosphere for different initial spin inclination is studied in the next section.
Fig. 10 Evolution of a SunVenuslike 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/s^{2}) 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 Venuslike 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 Venuslike 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 presentday 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 stableequilibrium 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 prograde 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 highinclination state and keep the rotation on the prograde 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 highinclination states. None of the paths shown in Fig. 13b crosses the null spin state. Any positive rotation with a lowinclination 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 Sunlike star over the time from the beginning to the end of the MS. The dashdotted 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 Venuslike planet. The luminosity evolution of the Sunlike star in ESPEM comes from simulations with the stellar evolution code STAREVOL (Amard et al. 2016). Figure 14 shows the luminosity variation of the Sunlike star we considered. The thermal Love number was computed with the luminosity dependence with the formulation of AuclairDesrotour et al. (2017b) as (13)
with L_{★} and M_{★} the stellar mass and luminosity respectively, R the radius of the planet, a the semimajor 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, ω_{0} the radiative frequency, T_{0} 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 nonevolving.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 .
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 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 Venuslike 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 semimajor 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_{2} associated with its atmosphere and internal structure are known.
First we focused in Sect. 3.1 on the eccentricitydriven 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 inclinationdriven 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 Venuslike 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 despinning 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 coremantle 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 Maxwelllike frequencydependent model of thermal tides overestimates the strength of the thermal tides at high frequencies, as the model was fit for the lowfrequency regime (AuclairDesrotour 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) threebody 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 eccentricitydriven SORs and also inclinationdriven 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 obliquitydriven 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 AuclairDesrotour 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, , ω, and ε the semimajor 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 spinaxis 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, , with r_{g} the gyration radius (which represents the internal density distribution). β is the reduced mass β = M_{p}M_{★}/(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 (with and n the spin rate and the mean motion respectively). The phases v_{lmpq} are defined as , with ℳ the mean anomaly. k_{l} is the modulus of the complex Love number of degree l. The tidal potential energy V_{1} is defined from the Hamiltonian formalism as (Boué & Efroimsky 2019) ℋ = ℋ + V_{1}. The tidal perturbing potential ℛ is computed with ℛ = V_{1}/β. The perturbing potential U (Eq. 1) within the DarwinKaula formalism can be related to R with (Boué & Efroimsky 2019). Then, the tidal perturbing potential ℛ is expressed in the formalism of Kaula (1964) as (A.1)
where the subscript ′ (i.e., r′, a′, e′, i′, , and ω′) corresponds to the coordinate of the perturber, and the parameters without subscript (i.e. r, a, e, i, , and ω) correspond to the coordinate where the tides are evaluated. In the following, both the perturber and the body on which the tides retroact 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) (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 (A.3)
The eccentricity derivative equation is computed from the Hamiltonian equation as (A.4)
Then, the equation of the eccentricity derivative is computed with the eccentricity squared e^{2} to avoid singularities when the eccentricity tends to zero. Thus, we find (A.5)
The equation of the inclination derivative is also defined from the Hamiltonian formalism as (A.6)
Then, we carry on the derivative of the perturbing function ℛ, and the equation can be written as (A.7)
The equation of the longitude of the ascending node derivative is given as (A.8)
Then, after some algebra, the equation can be written as (A.9)
The equation of the argument of the periastron derivative is defined as (A.10)
Then, the equation can be written as (A.11)
Finally, the equation of the inclination of the spin axis is defined as (A.12)
Then, after some algebra, we find (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 (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.
Inclination polynomials F_{lmp}(i) for l = 2 (Kaula 1964)
Eccentricity polynomials G_{lpq}(e) from Cayley (1861), up to order 7 in eccentricity.
References
 Ahuir, J., Strugarek, A., Brun, A. S., & Mathis, S. 2021, A & A, 650, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Amard, L., Palacios, A., Charbonnel, C., Gallet, F., & Bouvier, J. 2016, A & A, 587, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Andrade, E. N. D. C. 1910, Proc. Roy. Soc. Lond. Ser. A, 84, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Andrade, E. N. D. C. 1914, Proc. Roy. Soc. Lond. Ser. A, 90, 329 [NASA ADS] [CrossRef] [Google Scholar]
 Armann, M., & Tackley, P. J. 2012, J. Geophys. Res. (Planets), 117, E12003 [Google Scholar]
 AuclairDesrotour, P., Laskar, J., & Mathis, S. 2017a, A & A, 603, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 AuclairDesrotour, P., Laskar, J., Mathis, S., & Correia, A. C. M. 2017b, A & A, 603, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 AuclairDesrotour, P., Leconte, J., & Mergny, C. 2019, A & A, 624, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bagheri, A., Efroimsky, M., CastilloRogez, J., et al. 2022, Adv. Geophys., 63, 231 [NASA ADS] [CrossRef] [Google Scholar]
 Benbakoura, M., Réville, V., Brun, A. S., Le PoncinLafitte, C., & Mathis, S. 2019, A & A, 621, A124 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bills, B. G., Navarro, T., Schubert, G., Ermakov, A., & Górski, K. M. 2020, Icarus, 340, 113568 [NASA ADS] [CrossRef] [Google Scholar]
 Bolmont, E., & Mathis, S. 2016, Celest Mech Dyn Astr, 126, 275 [NASA ADS] [CrossRef] [Google Scholar]
 Bolmont, E., Breton, S. N., Tobie, G., et al. 2020a, A & A, 644, A165 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bolmont, E., Demory, B. O., BlancoCuaresma, S., et al. 2020b, A & A, 635, A117 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Boué, G., Correia, A. C. M., & Laskar, J. 2016, Celest. Mech. Dyn. Astron., 126, 31 [CrossRef] [Google Scholar]
 Boué, G., & Efroimsky, M. 2019, Celest. Mech. Dyn. Astron., 131, 30 [CrossRef] [Google Scholar]
 Bower, D. J., Kitzmann, D., Wolf, A. S., et al. 2019, A & A, 631, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 CastilloRogez, J. C., Efroimsky, M., & Lainey, V. 2011, J. Geophys. Res. (Planets), 116, E09008 [CrossRef] [Google Scholar]
 Cayley, A. 1861, MmRAS, 29, 191 [NASA ADS] [Google Scholar]
 Chambers, J. E. 2004, Earth Planet. Sci. Lett., 223, 241 [CrossRef] [Google Scholar]
 Chapman, S., & Lindzen, R. 1970, Atmospheric Tides. Thermal and Gravitational [Google Scholar]
 Connolly, J. A. D. 2005, Earth Planet. Sci. Lett., 236, 524 [CrossRef] [Google Scholar]
 Correia, A. C. M., & Laskar, J. 2001, Nature, 411, 767 [Google Scholar]
 Correia, A. C. M., & Laskar, J. 2003a, J. Geophys. Res. (Planets), 108, 5123 [NASA ADS] [CrossRef] [Google Scholar]
 Correia, A. C. M., & Laskar, J. 2003b, Icarus, 163, 24 [CrossRef] [Google Scholar]
 Correia, A. C., Laskar, J., & de Surgy, O. N. 2003, Icarus, 163, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Correia, A. C. M., Boué, G., Laskar, J., & Rodríguez, A. 2014, A & A, 571, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Darwin, G. H. 1879, The Observatory, 3, 79 [NASA ADS] [Google Scholar]
 Dobrovolskis, A. R., & Ingersoll, A. P. 1980, Icarus, 41, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Dumoulin, C., Tobie, G., Verhoeven, O., Rosenblatt, P., & Rambaux, N. 2017, J. Geophys. Res. (Planets), 122, 1338 [NASA ADS] [CrossRef] [Google Scholar]
 Edwards, B., & Tinetti, G. 2022, AJ, 164, 15 [NASA ADS] [CrossRef] [Google Scholar]
 Efroimsky, M. 2012a, Celest. Mech. Dyn. Astron., 112, 283 [NASA ADS] [CrossRef] [Google Scholar]
 Efroimsky, M. 2012b, ApJ, 746, 150 [Google Scholar]
 Efroimsky, M., & Makarov, V. V. 2013, ApJ, 764, 26 [NASA ADS] [CrossRef] [Google Scholar]
 Forget, F., Wordsworth, R., Millour, E., et al. 2013, Icarus, 222, 81 [NASA ADS] [CrossRef] [Google Scholar]
 Garvin, J. B., Getty, S. A., Arney, G. N., et al. 2022, Planet. Sci. J., 3, 117 [NASA ADS] [CrossRef] [Google Scholar]
 Gevorgyan, Y., Boué, G., Ragazzo, C., Ruiz, L. S., & Correia, A. C. 2020, Icarus, 343, 113610 [NASA ADS] [CrossRef] [Google Scholar]
 Gold, T., & Soter, S. 1969, Icarus, 11, 356 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P. 1966, AJ, 71, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & Soter, S. 1966, Icarus, 5, 375 [Google Scholar]
 Gomes, G. O., Bolmont, E., & BlancoCuaresma, S. 2021, A & A, 651, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Greene, T. P., Line, M. R., Montero, C., et al. 2016, ApJ, 817, 17 [NASA ADS] [CrossRef] [Google Scholar]
 Henning, W. G., O'Connell, R. J., & Sasselov, D. D. 2009, ApJ, 707, 1000 [NASA ADS] [CrossRef] [Google Scholar]
 Hut, P. 1981, A & A, 99, 126 [NASA ADS] [Google Scholar]
 Ingersoll, A. P., & Dobrovolskis, A. R. 1978, Nature, 275, 37 [CrossRef] [Google Scholar]
 Izsak, I. G., Gerard, J., Efimba, R., & Barnett, M. 1964, SAO Special Report, 140 [Google Scholar]
 Kaula, W. M. 1961, Geophys. J., 5, 104 [NASA ADS] [CrossRef] [Google Scholar]
 Kaula, W. M. 1964, Rev. Geophys. Space Phys., 2, 661 [Google Scholar]
 Laskar, J. 1990, Icarus, 88, 266 [NASA ADS] [CrossRef] [Google Scholar]
 Leconte, J., Forget, F., Charnay, B., Wordsworth, R., & Pottier, A. 2013a, Nature, 504, 268 [Google Scholar]
 Leconte, J., Forget, F., Charnay, B., et al. 2013b, A & A, 554, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Leconte, J., Wu, H., Menou, K., & Murray, N. 2015, Science, 347, 632 [NASA ADS] [CrossRef] [Google Scholar]
 Love, A. E. H. 1909, MNRAS, 69, 476 [NASA ADS] [Google Scholar]
 Makarov, V. V., & Efroimsky, M. 2013, ApJ, 764, 27 [Google Scholar]
 Makarov, V. V., Berghea, C. T., & Efroimsky, M. 2018, ApJ, 857, 142 [CrossRef] [Google Scholar]
 Mathis, S., & Le PoncinLafitte, C. 2009, A & A, 497, 889 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 McCarthy, C., & CastilloRogez, J. C. 2013, Astrophys. Space Sci. Lib., 356, 183 [NASA ADS] [CrossRef] [Google Scholar]
 Murray, C. D., & Dermott, S. F. 1999, Solar system Dynamics (Cambridge University Press) [Google Scholar]
 Ogilvie, G. I. 2014, ARA & A, 52, 171 [NASA ADS] [CrossRef] [Google Scholar]
 Remus, F., Mathis, S., Zahn, J. P., & Lainey, V. 2012, A & A, 541, A165 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Renaud, J. P., & Henning, W. G. 2018, ApJ, 857, 98 [Google Scholar]
 Smrekar, S. E., Hensley, S., Dyar, M. D., et al. 2020, in LPI Contribution, 2132, 51st Lunar and Planetary Science Conference, 1449 [Google Scholar]
 Tinetti, G., Eccleston, P., Haswell, C., et al. 2021, ArXiv eprints, [arXiv:2184.84824] [Google Scholar]
 Tisserand, F. 1889, Traité de mécanique céleste: Tome I, Perturbations des planètes d'après la méthode de la variation des constantes arbitraries (GauthierVillars) [Google Scholar]
 Tobie, G., Mocquet, A., & Sotin, C. 2005, Icarus, 177, 534 [Google Scholar]
 Tobie, G., Grasset, O., Dumoulin, C., & Mocquet, A. 2019, A & A, 630, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Turbet, M., Leconte, J., Selsis, F., et al. 2016, A & A, 596, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Walterová, M., & Behounková, M. 2020, ApJ, 900, 24 [CrossRef] [Google Scholar]
 Widemann, T., Ghail, R., Wilson, C. F., & Titov, D. V. 2020, in AGU Fall Meeting Abstracts, 2020, P02202 [Google Scholar]
 Wordsworth, R. 2015, ApJ, 806, 180 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Eccentricity polynomials G_{lpq}(e) from Cayley (1861), up to order 7 in eccentricity.
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 (substellar) 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 U_{lmpq} 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 dashdotted 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 as a function of the excitation frequency of a Venuslike 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 dashdotted 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 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 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 dashdotted 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 Venuslike 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 SunVenuslike 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/s^{2}) 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 Sunlike star over the time from the beginning to the end of the MS. The dashdotted 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.