11institutetext: Observatoire de Genève, Centre pour la Vie dans l’Univers,
Université de Genève, 51 Chemin Pegasis, 1290 Sauverny, Switzerland.
11email<EMAIL_ADDRESS><EMAIL_ADDRESS>22institutetext:
Laboratoire de Planétologie et Géosciences, UMR-CNRS 6112, Nantes Université,
2 rue de la Houssinière, BP 92208, 44322 Nantes Cedex 3, France
33institutetext: AIM, CEA, CNRS, Université Paris-Saclay, Université Paris
Diderot, Sorbonne Paris Cité, 91191 Gif-sur-Yvette Cedex, France.

# Spin evolution of Venus-like planets subjected
to gravitational and thermal tides

A. Revol 11 E. Bolmont 11 G. Tobie 22 C. Dumoulin 22 Y. Musseau 22 S. Mathis
33 A. Strugarek 33 A.S. Brun 33

(Received 26 December 2022 / Accepted 28 February 2023)

###### Abstract

Context. The arrival of powerful instruments will provide valuable data for
the characterization of rocky exoplanets. Rocky planets are mostly found in
close-in orbits. They are therefore usually close to the circular-coplanar
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 spin-orbit
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 Sun-Venus-like system in terms of stellar parameters, orbital configuration
and planet size and mass. The Sun-Venus 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 Venus-like exoplanets.

Results.  The formalism of Kaula associated with an Andrade rheology allows
spin orbit resonances on pure rocky worlds. Similarly to the high-order spin-
orbit resonances induced by eccentricity, the spin obliquity allows the
excitation of high-frequency Fourier modes that allow some spin-orbit
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
anti-synchronization 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 Venus-like 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: terretrial planets – planets and satellites: dynamical
evolution and stability – planet-star interactions

##  1 Introduction

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

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

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

In Section 2 we introduce the tidal model we used for the solid and thermal
tides and the implementation in the ESPEM code. In Section 3 we discuss the
evolution of the spin of a Venus-like planet when we only consider the
influence of the solid tide. In particular, we discuss the well-known
eccentricity-driven SORs and the less well-known inclination driven SORs. In
Section 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 Section 5.

##  2 New model of planetary tides in ESPEM

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

(a) Schematic representation of the atmospheric redistribution caused by the
stellar heating on a synchronously rotating planet. The arrows show the
movement of the particles of the atmosphere pushed from the hot spot (sub-
stellar) toward the cold spots (morning and evening spots).

(b) Delayed deformation of the atmosphere with respect to the position of the
star (angle $\delta_{a}$ in the schema) due to the rotation of the planet.
$\Omega$ is the spin rate of the planet, and is $n$ the mean motion and the
star.

(c) Two tidal contributions, gravitational and thermal, acting in opposition
on each other.

Figure 1: Tidal elongations of a rotating planet, composed of a solid core and
a gaseous atmosphere, submitted to both gravitational and thermal forcing.
Figures inspired from Correia & Laskar (2003).

In the following, we introduce in Section 2.1 the formalism we used (Kaula)
and in Section 2.2 the complex Love number, which allows us to express the
tidal potential of the deformed planet in terms of Fourier series. In Section
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 Section 2.4 we design a homogeneous interior model for the planet that
counterbalances the thermal tides for the excitation of Venus. In Section 2.5
we discuss the corresponding orbital and rotational equations to finally
described the implementation in the ESPEM code (Section 2.6).

###  2.1 Kaula formalism

In order to compute the tidal response of a body to a tidal perturbation, we
need to use a formalism that is general enough to encapsulate the frequency-
dependent response of a body. This response either requires a decomposition of
the tidal potential created by the perturber (hereafter perturbing potential)
into Fourier harmonic modes as developed by Kaula (1964), or a time-domain
approach as proposed by Correia et al. (2014) and Gevorgyan et al. (2020).
Both models allow a study of more complex and realistic rocky and icy bodies
(Efroimsky & Makarov, 2013; Bolmont et al., 2020a). We used the formalism
developed by Darwin (1879) and adapted by Kaula (1964), hereafter the Darwin-
Kaula formalism. The Darwin-Kaula theory of bodily tides provides the
expression of the perturbing potential of a disturbed body in Fourier series
as

$\begin{split}U=&\sum_{l=2}^{\infty}\sum_{m=0}^{l}\sum_{p=0}^{l}\sum_{q\in\mathbb{Z}}U_{lmpq}(a,e,i,\sigma_{lmpq}),\\\
\end{split}$ (1)

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

Figure 2: Schematic representation of the contribution of the tidal modes
$l,m,p,q=(2200),(2211),\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\
(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.)

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 Sec. 2.3). When the perturber is far enough away, we can keep the
development of the gravitational potential at the quadrupolar order only,
$l=2$ (Makarov & Efroimsky, 2013; Mathis & Le Poncin-Lafitte, 2009). We
restricted the eccentricity expansion numerically up to the order 7, which
corresponds to the index 7 in the summation over $q$ and eccentricities up to
$0.3$.

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

$U_{2mpq}^{\text{tidal
bulge}}(\sigma_{2mpq})=\bar{k_{2}}(\sigma_{2mpq})\leavevmode\nobreak\
U_{2mqp}^{\text{tidal perturbation}}(\sigma_{2mpq}),$ (2)

where $U_{2mpq}^{\text{tidal bulge}}(\sigma_{2mpq})$ corresponds to the
quadrupolar component of the potential Eq. 1 ($l=2$), and
$\bar{k_{2}}(\sigma_{2mpq})$ (hereafter $\bar{k_{2}}$) is related to the
amplitude of the complex quadrupolar Love number (Efroimsky, 2012a) at the
quadrupolar mode $l=2$ with a frequency dependence with $\sigma_{2mpq}$.

Physically, the quadrupolar Love number quantifies the response of a body
submitted to a periodical external perturbation of frequency $\sigma_{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

$\bar{k_{2}}=\Re\leavevmode\nobreak\ (\bar{k_{2}})+i\leavevmode\nobreak\
\Im\leavevmode\nobreak\ (\bar{k_{2}})\leavevmode\nobreak\
=|\bar{k_{2}}|\leavevmode\nobreak\ \text{exp}\big{(}-i\epsilon_{2}\big{)}.$
(3)

Thus, we can link the phase of the exponent $\epsilon_{2}$ with the angle
between the tidal bulge and the position of the perturber with
$\delta=\epsilon_{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)

$T^{\text{grav}}=\frac{3}{2}\frac{\mathcal{G}M_{\star}R_{p}^{5}}{a^{6}}\Im(\bar{k}_{2}),$
(4)

with $\mathcal{G}$ the gravitational constant, $M_{\star}$ the stellar mass,
$R_{p}$ the planetary radius, $a$ the semi-major axis and $\Im(\bar{k}_{2})$
the imaginary part of the Love number, which can be linked with the well-known
dissipation factor $Q_{2}$ and the Love number modulus $k_{2}$ with
$\Im(\bar{k}_{2}(\sigma))=-k_{2}(\sigma)/Q_{2}(\sigma)\text{Sign}(\sigma)$
(Goldreich & Soter, 1966; Ogilvie, 2014; Bagheri et al., 2022).

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

Figure 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.

The Love number can be computed with the method of Bolmont et al. (2020a).
Following Efroimsky (2012b), the complex tidal Love number
$\bar{k}^{\text{grav}}_{2}$ is given by

$\bar{k}^{\text{grav}}_{2}=\frac{3}{2}\frac{1}{1+A_{2}J/\bar{J}}\leavevmode\nobreak\
,$ (5)

with $J=1/\mu$ the unrelaxed compliance (with $\mu$ the unrelaxed elastic
shear modulus in Pa) and $A_{2}$ defined by

$A_{2}=\frac{57J^{-1}}{8\pi\mathcal{G}\rho^{2}R_{p}^{2}}\leavevmode\nobreak\
,$ (6)

with $\rho$ the density. $\bar{J}$ is the complex compliance of the material
and defined in the formalism of the Andrade rheology with (Castillo-Rogez et
al., 2011)

$\bar{J}=J+\beta(i\sigma)^{-\alpha}\Gamma(1+\alpha)-\frac{i}{\eta\sigma}\leavevmode\nobreak\
,$ (7)

with $\beta$ a factor that describes the intensity of anelastic friction in
the material, $\Gamma$ the Gamma function, $\eta$ the shear viscosity,
$\sigma$ the excitation frequency, and $\alpha$ an experimentally fit
parameter that which represents the frequency dependence of the transient
response. A value of $\alpha$ 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 Venus-like planet with different end-member
temperature profiles.Little is known about the interior structure of Venus.
This will likely improve with the upcoming ES EnVision mission (Widemann et
al., 2020) and the NASA DaVinci (Garvin et al., 2022) and VERITAS (Smrekar et
al., 2020) missions. In the meantime, we considered four possible structures.
We used one multilayer profile (referred to as the reference profile) with
Earth-like viscosity values as a reference, two other profiles with viscosity
values divided by 10 or multiplied by 100 relative to the reference profile
and one homogeneous profile (see Sec. 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 section 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
code222http://www.perplex.ethz.ch (Connolly, 2005), which uses a temperature
profile from Armann & Tackley (2012) together with the shear modulus profile
from the compositional model V1 of Dumoulin et al. (2017). The viscosity was
computed as a function of the temperature and pressure profiles as (Dumoulin
et al., 2017)

$\eta=\frac{1}{2}A_{0}^{-1}d^{2.5}\mathrm{exp}\Big{(}\frac{E_{a}+PV_{a}}{RT}\Big{)},$
(8)

with $E_{a}$ and $V_{a}$ the activation and volume energy, respectively, and
$A_{0}$ the pre-exponential factor, which are parameters from the Arrhenius
equation and depend on the material and $d$ the grain size. The parameters of
the dry olivine considered in the upper mantle are
$E_{a}=300\leavevmode\nobreak\ \mathrm{kJ}\leavevmode\nobreak\
\mathrm{mol}^{-1}$, $A_{0}=6.08\times 10^{-19}\leavevmode\nobreak\
\mathrm{Pa}^{-1}\leavevmode\nobreak\ \mathrm{s}^{-1}$, with a grain size
$d=0.68\leavevmode\nobreak\ \mathrm{mm}$.

Figure 4 shows the internal profiles of the shear modulus $\mu$ and the
viscosity $\eta$ for the homogeneous model (see Sec 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).

Figure 4:  Shear modulus and viscosity profiles for the multilayer reference
model Vref and homogeneous model considered here, shown as solid red and black
lines, respectively. The dotted and dash-dotted red lines represent the two
profiles V0.1 and V100 derived from the multilayer reference model with a
viscosity multiplied by x0.1 and x100 times, respectively. The viscosity
$\eta$ is computed as in Dumoulin et al. (2017) using Eq. 8 of this work.

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

Figure 5:  Imaginary part of the gravitational Love number
$\Im(k_{2}^{\text{grav}})$ as a function of the excitation frequency of a
Venus-like planet for different viscosity profiles. The multilayer reference
profile Vref derived from Armann & Tackley (2012) is shown as the solid red
line. The dotted and dash-dotted red lines represent the V0.1 and V100
profiles derived from the multilayer reference one presented in Fig 4. In blue
we present the imaginary Love number $\Im(k_{2}^{\text{thermal}})$ 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 Sec 2.4. The
vertical dotted black line represents the absolute value of the current
frequency state of Venus.

###  2.3 Thermal Love number

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

To calculate the complex moment of the surface pressure field $\delta
p_{s}^{2}$, we used the work of Leconte et al. (2015), who assumed a Maxwell-
like frequency dependence (Ingersoll & Dobrovolskis, 1978; Gold & Soter, 1969;
Auclair-Desrotour et al., 2017b). More realistic frequency dependences have
been proposed by Auclair-Desrotour et al. (2019) as a generic formulation and
a scaling law, adapted for $N_{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., 2013b; Forget et al.,
2013; Leconte et al., 2013a) specifically tuned for the case of Venus to
reproduce the amplitude of the thermal tides on Venus today and used this
point to fit an analytical Maxwell-like solution. The analytical formulation
of the module of the complex moment of the pressure field is expressed as

$\delta{p_{s}^{2}}=-\frac{q_{0}}{1+i\frac{\sigma}{2\omega_{0}}},$ (9)

such that the imaginary part can be written as

$\Im(\delta{p_{s}^{2}})=\frac{q_{0}\frac{\sigma}{2\omega_{0}}}{1+\left(\frac{\sigma}{2\omega_{0}}\right)^{2}},$
(10)

with $q_{0}$ the amplitude of the quadrupole term of the pressure field at
zero frequency, $\omega_{0}$ the radiative frequency, and $\sigma$ 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 $\omega_{0}=3.77\times 10^{-7}$
s-1 (Leconte et al., 2015). Figure 6 shows the amplitude of the pressure bulge
$|\delta{p_{s}^{2}}|$ 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.

Figure 6:  Amplitude of the pressure bulge $|\delta{p_{s}^{2}}|$ as a function
of the normalized forcing frequency $(\Omega-n)/\Omega$. 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.

The thermal Love number $\Im(k_{2}^{\text{thermal}})$ can be determined from
the complex moment of the surface pressure field $\Im\big{(}\delta
p_{s}^{2}\big{)}$ 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,
2003; Leconte et al., 2015; Auclair-Desrotour et al., 2017a)

$T^{\text{thermal}}=\sqrt{\frac{24\pi}{5}}\frac{M_{\star}}{M_{p}}\frac{R_{p}^{6}}{a^{3}}\leavevmode\nobreak\
\Im\big{(}\delta p_{s}^{2}\big{)}\leavevmode\nobreak\ ,$ (11)

with $M_{\star}$ the stellar mass, $M_{p}$ and $R_{p}$ the planetary mass and
radius respectively, and $a$ the semi-major axis. The thermal equivalent Love
number $\Im(k_{2})$ can be written by identification with the expression of
the solid torque (Eq. 4) as

$\Im(k_{2}^{\text{thermal}})=-\sqrt{\frac{32\pi}{15}}\frac{a^{3}R_{p}}{\mathcal{G}M_{\star}M_{p}}\Im\big{(}\delta
p_{s}^{2}\big{)}\leavevmode\nobreak\ ,$ (12)

We note that in contrast to $\Im(k_{2}^{\text{grav}})$, the thermal Love
number $\Im(k_{2}^{\text{thermal}})$ depends on the semi-major axis of the
planet as the intrinsic response of the atmosphere depends on the flux
received by the planet. Then, the two torques, gravitational and thermal, do
not have the same dependence as to the semi-major axis, which allows an
equilibrium point at which the two tides compensate for each other. The
dependence of the mass of the atmosphere is contained in the surface pressure
term. We highlight that a more massive atmosphere does not necessarily lead to
stronger atmospheric tides. For a more massive atmosphere, the atmospheric
layers are more opaque to the stellar flux. Thus, as less stellar flux reaches
the surface, the thermal tides are damped for a more massive atmosphere. This
effect strongly depends on the composition of the atmosphere and requires a
better model to be taken into account. A more massive atmosphere is not
investigated in this study. The imaginary Love number
$\Im(k_{2}^{\text{thermal}})$ 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.

###  2.4 Equilibrium state between gravitational and thermal tides

An equilibrium state between the gravitational and thermal tides can be
determined by comparing their imaginary Love numbers. The two tides compensate
for each other when the addition of the two imaginary Love numbers is $0$
(when the two absolute values are equal; see Fig. 5). Figure 5 shows that the
$\Im(k_{2}^{\text{grav}})$ corresponding to the multilayer reference profile
(solid red curve) is always higher in amplitude than the
$\Im(k_{2}^{\text{thermal}})$ (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
$\Im(k_{2}^{\text{thermal}})$ (solid blue curve) has a broad range of
frequencies (from $10^{-7}s^{-1}$ to $2\times 10^{-5}\text{s}^{-1}$), where it
dominates the gravitational Love number $\Im(k_{2}^{\text{grav}})$ 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}\leavevmode\nobreak\ $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\times 10^{-5}\leavevmode\nobreak\ $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 high-
frequency regime. In our case, the Maxwell-like frequency-dependent model of
thermal tides overestimates the strength of the thermal tides at high
frequencies, as the model was fit for the low-frequency regime (Auclair-
Desrotour et al., 2019). We therefore consider our approach to be valid in the
low-frequency range, that is, for $|\Omega/n|<15$. As this equilibrium point
is very far from the present-day rotation of Venus, it would require a more
complex model that is beyond the scope of this study, and we did not
investigate it further. A better model, such as the parameterized model
proposed by Auclair-Desrotour et al. (2019), will be studied in the future.

Figure 7:  Spin derivative as a function of the rotation (in terms of
$\Omega/n$, $\Omega$ and $n$ the planetary spin and mean motion respectively).
In the top panel, the red lines correspond to the solid tides, and the blue
lines correspond to the cases with solid and atmospheric tides. The solid
lines correspond to the reference multilayer profile Vref. The dotted and
dash-dotted lines correspond to the V0.1 and V100 profiles, respectively. In
the bottom panel, the green line represents the solid tides associated with
the homogeneous body (see Sec. 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).

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 Sec 2.1 (Eq. 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 $\alpha=0.25$
and only fit the value of the viscosity parameter $\text{log}(\eta)$. The
other parameters are the rigidity $\text{log}(\mu)=10.02$ in Pa and the ratio
of the Andrade and Maxwell time $\tau_{A}/\tau_{M}=0.89$ (Castillo-Rogez et
al., 2011). The homogeneous profile that fits the thermal Love number
$\Im(k_{2}^{\text{thermal}})$ at the Venus frequency is found with
$\text{log}(\eta)=22.18$ and is plotted in Fig. 5 (green curve). Fig 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.

###  2.5 Secular equations

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

In this formalism, the inclination is defined by the angle between the orbital
plane and the equatorial plane, which in other words corresponds to the angle
between the orbital angular momentum and the planetary spin angular
momentum333This angle is also often referred to as the obliquity, for
instance, the obliquity of the Earth is about 23 degrees. In this study, we
would thus say that the inclination of the Earth is 23 degrees..

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

###  2.6 Implementation in the ESPEM code

We implemented the secular equations of Boué & Efroimsky (2019) in the code
ESPEM (Benbakoura et al., 2019; Ahuir et al., 2021). This is a secular code
integrating the dynamical evolution of a star-planet system. The code takes
the coupling between the two layers of low-mass stars into account (convective
and radiative layers), as well as the effect of the stellar wind and the
torque due to the tides raised by the planet on the 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 semi-major axis, we also implemented the equations governing
further osculating elements of the planet, such as the eccentricity, orbital
inclination, longitude of ascending node, argument of periapsis, and planetary
spin. The equations for spin, eccentricity, and inclination, longitude of
ascending node and argument of periastron can be found in the Appendix A, Eq.
15 to 19.

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 Sec 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 Sec 2.4. As explained in Sec 2.1, several
frequencies are excited depending on the eccentricity and inclination. These
frequencies were computed for each time step. The real and imaginary parts of
the Love number were interpolated linearly from their frequency dependence in
the data file. These interpolated values were then used to compute the
derivatives of all the quantities mentioned before.

We used the parameters of a Sun-Venus system, that is, a Venus-like mass and
radius planet orbiting at $0.723$ AU. The parameters we used are listed in
table 1. The initial spin rate is shown from $\Omega/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.

Table 1: Numerical values used in the case of a Sun-Venus-like system. Parameter | Values
---|---
Star Mass ($M_{\odot}$) | 1
Planet Mass ($M_{\text{Earth}}$) | $0.815$
Planet Radius ($R_{\text{Earth}}$) | $0.857$
Semi-major axis (AU) | $0.723$
Eccentricity | $\\{0,0.1,0.2\\}$
Spin Inclination (degrees) | $\\{0,5,50,120,130\\}$

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}$ Gyrs order of magnitude(e.g., Bolmont & Mathis,
2016).

##  3 Impact of the gravitational tides alone

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

###  3.1 Eccentricity-driven spin-orbit resonances

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

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

The order of the resonance in which the planet capture is shown when we plot
the spin derivative as a function of $\Omega/n$ ($\Omega$ and $n$ the
planetary spin and mean motion respectively). Figure 9 shows how this quantity
evolves with $\Omega/n$ for a fixed eccentricity (Fig. 9a) and how it evolves
with $\Omega/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\Omega/dt=0$, and therefore, only one possible equilibrium
for the rotation state. The equilibrium is centered at $\Omega/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 $\Omega/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 higher-order SORs. The
synchronization is accessible with $e=0$, while the 3:2 SOR becomes accessible
at $e=0.06$, the 2:1 SOR at $0.16$, and the 5:2 SOR at $0.26$. The eccentric
cases of the Fig. 9a are overplotted in the color map with the two horizontal
dotted black lines. The evolutions shown in Fig. 8 are plotted in Fig. 9b and
9c with the three colored arrows (with identical colors in the two figures).
The spin quickly decreases in the SOR associated with its eccentricity.
Because the eccentricity is damped by the tides, the SORs remain until the
eccentricity becomes too low to stably maintain these configurations. Figure
9c represents the eccentricity derivative map in the eccentricity versus
rotation state plane. In red we show the area in which the eccentricity
increases, and in blue the area in which it decreases. In particular, the
eccentricity appears to be slightly excited for the $e=0.1$ case shown in Fig.
8. This behavior can be explained with the derivative map of Fig. 9c. The
eccentricity in the $e=0.2$ case of Fig. 8 also appears to be excited before
the state when the rotation reached the $2:1$ SOR.

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

Figure 8:  Evolution of a Venus-like planet with an initial rotation of $100$
days for three initial eccentricities (i.e., the null eccentricity, and the
$0.1$ and $0.2$ eccentricities). The top panel shows the evolution of the
planet rotation rate in terms of $\Omega/n$ ($\Omega$ and $n$ are the spin and
mean motion, respectively). The bottom panel shows the variation in
eccentricity $\Delta e$ (i.e., $(e-e(0))/e(0)$).  Figure 9:  Spin derivative
$\mathrm{d}\Omega/\mathrm{d}t$ ($rad/s^{2}$). The left panel shows the spin
derivative as a function of the rotation state $\Omega/n$ ($\Omega$ 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 $\Omega/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 $\mathrm{d}\Omega/\mathrm{d}t=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.

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

###  3.2 Inclination-driven spin-orbit resonances

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

Figure 10: Evolution of a Sun-Venus-like system. The top panel shows the
evolution of the rotation state in terms of $\Omega/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.

Figure 10 shows the evolution of a Venus-like planet with the multilayer
internal reference structure (see Sec. 2.2) after $2\times 10^{8}$ years 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\times 10^{7}$ years, 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 $\Omega/n$ and the inclination without eccentricity.
Figure 11b shows the spin derivative strength in the plane inclination versus
rotation state $\Omega/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 $\Omega/n$, and inversely for negative values of
$\Omega/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 $\Omega/n=-2$. A symmetry
with respect to a $90$-degree inclination exists. This symmetry is clearly
visible in the middle panel of Fig 11b. In particular, the torque at
$130$-degree inclination is symmetric of the $50$-degree inclination. We
highlight that no SORs lie above the 2:1 SOR in rotation. Figure 11b shows no
SORs close to the $\Omega/n=3\leavevmode\nobreak\
\text{or}\leavevmode\nobreak\ =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 Fig 11b-11c). This resonance island is stable for inclinations
greater than $15$ and lower than $120$ degrees. Thus, if the initial spin is
higher than $\Omega/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 $\mathrm{d}i/\mathrm{d}t$ plotted in Fig. 12. It shows a positive-
inclination derivative for a spin higher than the 2:1 SOR and an inclination
lower than about $80$ degrees (bottom right corner of the figure). It also
shows that the inclination should also increase when the spin is slightly
higher than the synchronization and for an inclination lower than 100 degrees
(red area in the vicinity of the 1:1 SOR). This behavior is absent in Fig. 10
because the spin reaches the synchronization very quickly.

Figure 11:  Spin derivative $\mathrm{d}\Omega/\mathrm{d}t$ (in $rad/s^{2}$).
The left panel shows the spin derivative as a function of the rotation state
$\Omega/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 $\Omega/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 Sec 2.5).

Higher inclination cases can show an interesting behavior. Fig. 11b clearly
shows that for a high initial inclination of about $105$ degrees, the
synchronization state (i.e, $\Omega/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
$\Omega/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
$\Omega/n=0$ axis and on the $90$-degree axis of the derivative map (Fig.
11b).

Figure 12:  Same as Fig 11b, with the inclination derivative
$\mathrm{d}i/\mathrm{d}t=f(\Omega/n,i)$, in $rad/s$, instead of the spin
derivative $\mathrm{d}\Omega/\mathrm{d}t$.

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

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

##  4 Venus-like atmospheric tides

As previous studies showed, the current spin state of Venus cannot be
reproduced by involving the solid tides alone (Gold & Soter, 1969;
Dobrovolskis & Ingersoll, 1980; Correia & Laskar, 2001, 2003; Correia et al.,
2003; Correia & Laskar, 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$ hr, 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 (Sec 4.1) explores the evolution of a Venus-like planet in
the spin and inclination parameter space, with a nonevolving atmosphere, a
constant luminosity, and a nonevolving internal profile. Section 4.2 explores
the effect of the luminosity evolution of the host star, accounting for a
simple prescription for the atmospheric evolution.

###  4.1 Constant luminosity, nonevolving atmosphere

Leconte et al. (2015) fit the parameters of their analytical solution of the
pressure bulge (Eq. 9) to their GCM simulation to model the thermal tides.
These parameters are given in Section 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 Section 2.2) and the fitted homogeneous model
(introduced in Section 2.4). For the thermal tides, we used the analytical
model of thermal tides fit on the present-day Venus (introduced in Section
2.3.) The frequency dependence of the corresponding Love numbers is given in
Fig. 5.

As discussed in Section 2.4, the solid tides corresponding to the multilayer
model and its two variants (the V0.1 and V100 profiles) do not allow a stable-
equilibrium point close to the current frequency of Venus as they are either
too strong or too weak. We fit a homogeneous hot profile, using the method of
Bolmont et al. (2020a), in order to find an equilibrium point close to the
Venusian frequency (see Sec 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 Sec 3.2, we used the derivative maps of $d\Omega/dt$ and $di/dt$ as a
function of $\Omega/n$ and $i$ to represent the evolution of the system in Fig
13. As discussed in Sec 2.4, we constrained our study at low spin rates.
Because no inclined SORs are higher than the 2:1 (see Sec 3.2), the initial
spin rate of the simulations was set to $\Omega/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). Figure 13a and 13b 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 high-inclination 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, 2003).

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

Figure 13:  Left panel: Inclination derivative $di/dt$ (in $rad/s$) as a
function of the spin $\Omega/n$ ($\Omega$ and $n$ are the planetary spin and
mean motion respectively) and for different inclinations (from $0$ to $180$
degrees). Right panel: Spin derivative $d\Omega/dt$ (in $rad/s^{2}$) as a
function of the spin $\Omega/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 (left panel) and the spin derivative (right panel). The
dotted orange lines represent the null points of the inclination derivative
from the left panel (in red in the left panel). 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.2 Luminosity variation

Dynamical studies of the thermal tides (Correia & Laskar, 2001, 2003; Leconte
et al., 2015) have considered a constant luminosity. As the thermal forcing
depends on the heat flux of the host star, we also investigated the effect of
an evolving luminosity on the rotation evolution of a Venus-like planet. The
luminosity evolution of the Sun-like star in ESPEM comes from simulations with
the stellar evolution code STAREVOL (Amard et al., 2016). Figure 14 shows the
luminosity variation of the Sun-like star we considered. The thermal Love
number was computed with the luminosity dependence with the formulation of
Auclair-Desrotour et al. (2017b) as

$\begin{split}\Im(k_{2}^{atm}(\sigma))&=-\frac{4}{32}\frac{\kappa\tau\varsigma\epsilon
L_{\star}a}{R_{A}T_{0}M_{\star}R}\frac{\sigma}{\sigma^{2}+\omega^{2}_{0}},\\\
\end{split}$ (13)

with $L_{\star}$ and $M_{\star}$ the stellar mass and luminosity respectively,
$R$ the radius of the planet, $a$ the semi-major axis, $\tau$ a weight
parameter that gives the efficiency of the coupling between the atmosphere and
the surface ($0<\tau<1$), $\varsigma$ a shape factor depending on the spatial
distribution of tidal heat sources, $\kappa$ the power per mass unit radiated
by the atmosphere (where the atmosphere is assumed to behave like a graybody,
i.e., Newtonian cooling), $\epsilon$ the effective fraction of power absorbed
by the atmosphere, $\sigma$ the excitation frequency, $\omega_{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}/\mathcal{M}_{A}$ ($R_{GP}$ and $\mathcal{M}_{A}$ being the
perfect gas constant and the mean molar mass respectively), and $\alpha$ the
shape factor depending on the spatial distribution of tidal heat sources. The
values of the parameters we used are presented in table 2.

Table 2: Numerical values for the thermal Love number of Eq. LABEL:eq:thermal_Lovenumber_auclair. Parameter | Values | Units
---|---|---
$T_{0}$ | $737$ | K
$\tau$ | $1$ | -
$\varsigma$ | $0.19$ | -
$\kappa$ | $0.286$ | -
$\epsilon$ | $0.04$ | -
$\omega_{0}$ | $3.77E-7$ | s-1
$R_{GP}$ | $8.314$ | J mol-1 K-1
$\mathcal{M}_{A}$ | $43.45$ | g mol-1
Figure 14:  Luminosity variation of a Sun-like star over the time from the
beginning to the end of the MS. The dash-dotted line represents the simulation
start time. The dotted lines represent the current time state.

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

Figure 16 shows the spin versus inclination maps with the ESPEM simulation
overplotted, in the same manner as the Fig 13a and 13b. Figures 16a to 16h
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 panel 16b to panel 16h). 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 figure 11c. Figure
16d 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 16e to 16h
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 $\Omega_{\text{eq}}$, which satisfies
$k_{2}^{\text{grav}}(\sigma_{\text{eq}})=k_{2}^{\text{thermal}}(\sigma_{\text{eq}})$.

Figure 15 shows the evolution of the equilibrium rotation rate in terms of
$\Omega_{\text{eq}}/n$ and the evolution of the rotation rate $\Omega/n$ from
the ESPEM simulations shown in Fig. 16. The equilibrium states were determined
numerically by finding the rotation state $\Omega_{eq}$ that verifies the
equality between the gravitational and the thermal Love number
$\Im(k_{2}^{\text{grav}}(\Omega_{eq}/n))=\Im(k_{2}^{\text{thermal}}(\Omega_{eq}/n))$
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. 16. 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 16h shows that the thermal tides eventually
become stronger than the gravitational tides across a large parameter space as
the luminosity increases.

Figure 15:  Evolution of the rotational equilibrium state $\Omega_{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. 16.

Figure 16: 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 Sec. 4.2.

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).

##  5 Conclusion

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

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

In Sec 4.1 we investigated the effect of a thick Venus-like atmosphere with
the implementation of the analytical model of Leconte et al. (2015). We used
their fit parameters, chosen so that GCM simulations reproduce the current
state of Venus. We fit a homogeneous internal structure that allowed
gravitational tides to balance thermal tides at the frequency of Venus. We
must emphasize that the de-spinning of Venus is a difficult task. Then, we
constrained our work to lower initial rotation rates. We find that depending
on the initial spin rate and initial spin inclination, either a spin
inclination of about zero in the synchronization state or a state close to the
current retrograde rotation of Venus results, with a spin rate close to the
synchronization and a spin inclination of $180$ degrees. The synchronization
state (1:1 SOR) is reached when the planet starts with a spin inclination
lower than about $120$ degrees and a prograde spin in our simulations. The
latter state can be reached when the planet starts either with a high spin
inclination (higher than about $120$ degrees) and a prograde rotation, or with
a spin inclination lower than about $60$ degrees and a retrograde rotation.
Our results are consistent with the final spin state of Venus found by Correia
& Laskar (2001), who computed the evolution of the spin and obliquity of Venus
under the solar tides. We point out, however, that Correia & Laskar (2001)
used different models for the gravitational tides, which less appropriate than
an Andrade rheology for silicate bodies (i.e., the CTL, Hut, 1981; Goldreich,
1966; Efroimsky & Makarov, 2013) and they included the core-mantle friction.
The core-mantle friction helps to damp the spin inclination of the planet
(Correia & Laskar, 2001) and should help the gravitational tides to balance
the thermal tides. Furthermore, they also accounted for the chaotic motion in
the Solar System (Laskar, 1990). In particular, they reported that the chaotic
motion helps to transition from low to high inclination. We cannot reproduce
this with only two bodies. Further developments of the ESPEM code will include
these effects.

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

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

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

###### Acknowledgements.

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

## References

  * Ahuir et al. (2021) Ahuir, J., Strugarek, A., Brun, A. S., & Mathis, S. 2021, A&A, 650, A126
  * Amard et al. (2016) Amard, L., Palacios, A., Charbonnel, C., Gallet, F., & Bouvier, J. 2016, A&A, 587, A105
  * Andrade (1910) Andrade, E. N. D. C. 1910, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character, 84, 1
  * Andrade (1914) Andrade, E. N. D. C. 1914, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character, 90, 329
  * Armann & Tackley (2012) Armann, M. & Tackley, P. J. 2012, Journal of Geophysical Research (Planets), 117, E12003
  * Auclair-Desrotour et al. (2017a) Auclair-Desrotour, P., Laskar, J., & Mathis, S. 2017a, A&A, 603, A107
  * Auclair-Desrotour et al. (2017b) Auclair-Desrotour, P., Laskar, J., Mathis, S., & Correia, A. C. M. 2017b, A&A, 603, A108
  * Auclair-Desrotour et al. (2019) Auclair-Desrotour, P., Leconte, J., & Mergny, C. 2019, A&A, 624, A17
  * Bagheri et al. (2022) Bagheri, A., Efroimsky, M., Castillo-Rogez, J., et al. 2022, Advances in Geophysics, 63, 231
  * Benbakoura et al. (2019) Benbakoura, M., Réville, V., Brun, A. S., Le Poncin-Lafitte, C., & Mathis, S. 2019, A&A, 621, A124
  * Bills et al. (2020) Bills, B. G., Navarro, T., Schubert, G., Ermakov, A., & Górski, K. M. 2020, Icarus, 340, 113568
  * Bolmont et al. (2020a) Bolmont, E., Breton, S. N., Tobie, G., et al. 2020a, A&A, 644, A165
  * Bolmont et al. (2020b) Bolmont, E., Demory, B. O., Blanco-Cuaresma, S., et al. 2020b, A&A, 635, A117
  * Bolmont & Mathis (2016) Bolmont, E. & Mathis, S. 2016, 126, 275
  * Boué et al. (2016) Boué, G., Correia, A. C. M., & Laskar, J. 2016, Celestial Mechanics and Dynamical Astronomy, 126, 31
  * Boué & Efroimsky (2019) Boué, G. & Efroimsky, M. 2019, Celestial Mechanics and Dynamical Astronomy, 131, 30
  * Bower et al. (2019) Bower, D. J., Kitzmann, D., Wolf, A. S., et al. 2019, A&A, 631, A103
  * Castillo-Rogez et al. (2011) Castillo-Rogez, J. C., Efroimsky, M., & Lainey, V. 2011, Journal of Geophysical Research (Planets), 116, E09008
  * Cayley (1861) Cayley, A. 1861, MmRAS, 29, 191
  * Chambers (2004) Chambers, J. E. 2004, Earth and Planetary Science Letters, 223, 241
  * Chapman & Lindzen (1970) Chapman, S. & Lindzen, R. 1970, Atmospheric tides. Thermal and gravitational
  * Connolly (2005) Connolly, J. A. D. 2005, Earth and Planetary Science Letters, 236, 524
  * Correia et al. (2003) Correia, A. C., Laskar, J., & de Surgy, O. N. 2003, 163, 1
  * Correia et al. (2014) Correia, A. C. M., Boué, G., Laskar, J., & Rodríguez, A. 2014, A&A, 571, A50
  * Correia & Laskar (2001) Correia, A. C. M. & Laskar, J. 2001, Nature, 411, 767
  * Correia & Laskar (2003) Correia, A. C. M. & Laskar, J. 2003, Journal of Geophysical Research (Planets), 108, 5123
  * Correia & Laskar (2003) Correia, A. C. M. & Laskar, J. 2003, 163, 24
  * Darwin (1879) Darwin, G. H. 1879, The Observatory, 3, 79
  * Dobrovolskis & Ingersoll (1980) Dobrovolskis, A. R. & Ingersoll, A. P. 1980, Icarus, 41, 1
  * Dumoulin et al. (2017) Dumoulin, C., Tobie, G., Verhoeven, O., Rosenblatt, P., & Rambaux, N. 2017, Journal of Geophysical Research (Planets), 122, 1338
  * Edwards & Tinetti (2022) Edwards, B. & Tinetti, G. 2022, AJ, 164, 15
  * Efroimsky (2012a) Efroimsky, M. 2012a, Celestial Mechanics and Dynamical Astronomy, 112, 283
  * Efroimsky (2012b) Efroimsky, M. 2012b, ApJ, 746, 150
  * Efroimsky & Makarov (2013) Efroimsky, M. & Makarov, V. V. 2013, ApJ, 764, 26
  * Forget et al. (2013) Forget, F., Wordsworth, R., Millour, E., et al. 2013, Icarus, 222, 81
  * Garvin et al. (2022) Garvin, J. B., Getty, S. A., Arney, G. N., et al. 2022, The Planetary Science Journal, 3, 117
  * Gevorgyan et al. (2020) Gevorgyan, Y., Boué, G., Ragazzo, C., Ruiz, L. S., & Correia, A. C. 2020, Icarus, 343, 113610
  * Gold & Soter (1969) Gold, T. & Soter, S. 1969, Icarus, 11, 356
  * Goldreich (1966) Goldreich, P. 1966, AJ, 71, 1
  * Goldreich & Soter (1966) Goldreich, P. & Soter, S. 1966, Icarus, 5, 375
  * Gomes et al. (2021) Gomes, G. O., Bolmont, E., & Blanco-Cuaresma, S. 2021, A&A, 651, A23
  * Greene et al. (2016) Greene, T. P., Line, M. R., Montero, C., et al. 2016, ApJ, 817, 17
  * Henning et al. (2009) Henning, W. G., O’Connell, R. J., & Sasselov, D. D. 2009, ApJ, 707, 1000
  * Hut (1981) Hut, P. 1981, A&A, 99, 126
  * Ingersoll & Dobrovolskis (1978) Ingersoll, A. P. & Dobrovolskis, A. R. 1978, Nature, 275, 37
  * Izsak et al. (1964) Izsak, I. G., Gerard, J., Efimba, R., & Barnett, M. 1964, SAO Special Report, 140
  * Kaula (1961) Kaula, W. M. 1961, Geophysical Journal, 5, 104
  * Kaula (1964) Kaula, W. M. 1964, Reviews of Geophysics and Space Physics, 2, 661
  * Laskar (1990) Laskar, J. 1990, Icarus, 88, 266
  * Leconte et al. (2013a) Leconte, J., Forget, F., Charnay, B., Wordsworth, R., & Pottier, A. 2013a, Nature, 504, 268
  * Leconte et al. (2013b) Leconte, J., Forget, F., Charnay, B., et al. 2013b, A&A, 554, A69
  * Leconte et al. (2015) Leconte, J., Wu, H., Menou, K., & Murray, N. 2015, Science, 347, 632
  * Love (1909) Love, A. E. H. 1909, MNRAS, 69, 476
  * Makarov et al. (2018) Makarov, V. V., Berghea, C. T., & Efroimsky, M. 2018, ApJ, 857, 142
  * Makarov & Efroimsky (2013) Makarov, V. V. & Efroimsky, M. 2013, ApJ, 764, 27
  * Mathis & Le Poncin-Lafitte (2009) Mathis, S. & Le Poncin-Lafitte, C. 2009, A&A, 497, 889
  * McCarthy & Castillo-Rogez (2013) McCarthy, C. & Castillo-Rogez, J. C. 2013, in Astrophysics and Space Science Library, Vol. 356, Astrophysics and Space Science Library, ed. M. S. Gudipati & J. Castillo-Rogez, 183
  * Murray & Dermott (1999) Murray, C. D. & Dermott, S. F. 1999, Solar system dynamics (Cambridge university press)
  * Ogilvie (2014) Ogilvie, G. I. 2014, ARA&A, 52, 171
  * Remus et al. (2012) Remus, F., Mathis, S., Zahn, J. P., & Lainey, V. 2012, A&A, 541, A165
  * Renaud & Henning (2018) Renaud, J. P. & Henning, W. G. 2018, ApJ, 857, 98
  * Smrekar et al. (2020) Smrekar, S. E., Hensley, S., Dyar, M. D., et al. 2020, in LPI Contribution, Vol. 2132, 51st Lunar and Planetary Science Conference, 1449
  * Tinetti et al. (2021) Tinetti, G., Eccleston, P., Haswell, C., et al. 2021, arXiv e-prints, arXiv:2104.04824
  * Tisserand (1889) 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 arbitraires (Gauthier-Villars)
  * Tobie et al. (2019) Tobie, G., Grasset, O., Dumoulin, C., & Mocquet, A. 2019, A&A, 630, A70
  * Tobie et al. (2005) Tobie, G., Mocquet, A., & Sotin, C. 2005, Icarus, 177, 534
  * Turbet et al. (2016) Turbet, M., Leconte, J., Selsis, F., et al. 2016, A&A, 596, A112
  * Walterová & Běhounková (2020) Walterová, M. & Běhounková, M. 2020, ApJ, 900, 24
  * Widemann et al. (2020) Widemann, T., Ghail, R., Wilson, C. F., & Titov, D. V. 2020, in AGU Fall Meeting Abstracts, Vol. 2020, P022–02
  * Wordsworth (2015) Wordsworth, R. 2015, The Astrophysical Journal, 806, 180

##  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$, $\bar{\Omega}$, $\omega$, and $\varepsilon$ the semi-major
axis, the eccentricity, the inclination, the longitude of ascending node, the
argument of pericenter and the inclination of the spin axis respectively, and
$\dot{\theta}$ 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 $\varepsilon$ is defined as the inclination of the spin with respect
to the inertial frame.

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

Hereafter, the star is taken as the secondary (subscript ⋆), and the planet as
the primary (subscript p), $C$ is the inertia momentum, that is,
$C=M_{p}R_{p}^{2}r_{g}$, with $r_{g}$ the gyration radius (which represents
the internal density distribution). $\beta$ is the reduced mass
$\beta=M_{p}M_{\star}/(M_{p}+M_{\star})$ and $\mathcal{G}$ the universal
gravitational constant. $F_{lmp}(i)$ and $G_{lpq}(e)$ are the inclination and
eccentricity polynomials respectively (see Appendix B). The tidal frequency is
defined as $\sigma_{lmpq}=(l-2p+q)n-m\dot{\theta}$ (with $\dot{\theta}$ and
$n$ the spin rate and the mean motion respectively). The phases $\nu_{lmpq}$
are defined as $\nu_{lmpq}=(l-2p)\omega+(l-2p+q)\mathcal{M}+m\bar{\Omega}$,
with $\mathcal{M}$ the mean anomaly. $k_{l}$ is the modulus of the complex
Love number $|\bar{k_{l}}|$ of degree $l$. The tidal potential energy $V_{1}$
is defined from the Hamiltonian formalism as (Boué & Efroimsky 2019)
$\mathcal{H}=\mathcal{H_{0}}+V_{1}$. The tidal perturbing potential
$\mathcal{R}$ is computed with $\mathcal{R}=-V_{1}/\beta$. The perturbing
potential $U$ (Eq. 1) within the Darwin-Kaula formalism can be related to
$\mathcal{R}$ with $\mathcal{R}=-\frac{M_{\star}}{\beta}U$ (Boué & Efroimsky
2019). Then, the tidal perturbing potential $\mathcal{R}$ is expressed in the
formalism of Kaula (1964) as

$\begin{split}\mathcal{R}(\mathbf{r},\mathbf{r^{\prime}})&=\sum_{l=0}^{+\infty}\Big{(}\frac{R_{p}}{a^{\prime}}\Big{)}^{l+1}\frac{\mathcal{G}M^{\prime}}{a^{\prime}}\Big{(}\frac{R_{p}}{a^{\prime}}\Big{)}^{l}\sum_{m=0}^{l}\frac{(l-m)!}{(l+m)!}(2-\delta_{0,m})\\\
&\sum_{p=0}^{l}F_{lmp}(i^{\prime})\sum_{q=-\infty}^{\infty}G_{lpq}(e^{\prime})\sum_{h=0}^{l}F_{lmh}(i)\sum_{j=-\infty}^{\infty}G_{lhj}(e)\\\
&k_{l}(\sigma_{lmpq})\cos{\big{[}(\nu^{\prime}_{lmpq}-m\dot{\theta}^{\prime})-(\nu_{lmhj}-m\dot{\theta})-\varepsilon_{l}(\sigma_{lmpq})\big{]}},\end{split}$
(14)

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

The spin derivative equation is expressed from the Hamiltonian formalism as
(Boué & Efroimsky 2019)

$\begin{split}\frac{d^{2}\theta}{dt^{2}}\Bigg{|}_{l=2}=-\frac{\beta}{C}\frac{\partial\mathcal{R}}{\partial\bar{\Omega}}&\\\
\end{split}$ (15)

We carry out the derivative of the perturbing function $\mathcal{R}$, and the
secular equation of the spin derivative is written after some algebra as

$\begin{split}\frac{d^{2}\theta}{dt^{2}}\Bigg{|}_{l=2}=-\frac{\mathcal{G}M_{\star}^{2}R^{5}}{a^{6}C}&\sum_{m=0}^{2}m\frac{(2-m)!}{(2+m)!}(2-\delta_{0,m})\sum_{p=0}^{2}{F_{2mp}(i)}^{2}\\\
&\sum_{q=-7}^{7}{G_{2pq}(e)}^{2}\Im\big{(}k_{2}(\sigma_{2mpq})\big{)}\\\
\end{split}$ (16)

The eccentricity derivative equation is computed from the Hamiltonian equation
as

$\begin{split}\frac{de}{dt}\Bigg{|}_{l=2}&=\frac{1-e^{2}}{na^{2}e}\frac{\partial\mathcal{R}}{\partial\mathcal{M}}-\frac{\sqrt{1-e^{2}}}{na^{2}e}\frac{\partial\mathcal{R}}{\partial\omega}\\\
\end{split}$ (17)

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

$\begin{split}\frac{de^{2}}{dt}\Bigg{|}_{l=2}&=-2\sqrt{(1-e^{2})}\sqrt{G(M_{p}+M_{\star})}\frac{M_{\star}}{M_{p}}\Big{(}\frac{R_{p}^{5}}{a^{13/2}}\Big{)}\\\
&\sum_{m=0}^{2}\frac{(2-m)!}{(2+m)!}(2-\delta_{0,m})\sum_{p=0}^{2}F_{2mp}(i)^{2}\\\
&\sum_{q=-7}^{7}G_{2pq}(e)^{2}\Big{(}\sqrt{1-e^{2}}(2-2p+q)-(2-2p)\Big{)}\Im\big{(}k_{2}(\sigma_{2mpq})\big{)}\\\
\end{split}$ (18)

The equation of the inclination derivative is also defined from the
Hamiltonian formalism as

$\begin{split}\frac{di}{dt}\Bigg{|}_{l=2}&=\frac{\beta}{C\dot{\theta}\sin{i}}\Big{(}\frac{\partial\mathcal{R}}{\partial\omega}-\cos{i}\frac{\partial\mathcal{R}}{\partial\bar{\Omega}}\Big{)}\\\
&-\frac{1}{na^{2}\sqrt{1-e^{2}}\sin{i}}\Big{(}\frac{\partial\mathcal{R}}{\partial\bar{\Omega}}-\cos{i}\frac{\partial\mathcal{R}}{\partial\omega}\Big{)}\\\
\end{split}$ (19)

Then, we carry on the derivative of the perturbing function $\mathcal{R}$, and
the equation can be written as

$\begin{split}\frac{di}{dt}\Bigg{|}_{l=2}&=\frac{1}{\sin{i}}\frac{M_{\star}}{M_{p}}\Big{(}\frac{R_{p}}{a}\Big{)}^{5}\sum_{m=0}^{2}\frac{(2-m)!}{(2+m)!}(2-\delta_{0m})\times\\\
&\sum_{p=0}^{2}\Big{[}\frac{\beta
n^{2}a^{2}}{C\dot{\theta}}\Big{(}m\cos{i}-(2-2p)\Big{)}-\frac{n}{\sqrt{1-e^{2}}}\Big{(}(2-2p)\cos{i}-m\Big{)}\Big{]}\times\\\
&{F_{2mp}(i)}^{2}\sum_{q=-7}^{7}G_{2pq}(e)^{2}\Im\big{(}k_{2}(\sigma_{2mpq})\big{)}\\\
\end{split}$ (20)

The equation of the longitude of the ascending node derivative is given as

$\begin{split}\frac{d\bar{\Omega}}{dt}\Bigg{|}_{l=2}&=\Big{(}\frac{\beta\cos{i}}{C\dot{\theta}\sin{i}}-\frac{\beta\cos{\varepsilon}\cos{\bar{\Omega}}}{C\dot{\theta}\sin{\varepsilon}}+\frac{1}{na^{2}\sqrt{1-e^{2}}\sin{i}}\Big{)}\frac{\partial\mathcal{R}}{\partial
i}\\\
&+\frac{\beta\cos{\varepsilon}\sin{\bar{\Omega}}\cot{i}}{C\dot{\theta}\sin{\varepsilon}}\frac{\partial\mathcal{R}}{\partial\bar{\Omega}}-\frac{\beta\cos{\varepsilon}}{C\dot{\theta}\sin{\varepsilon}}\frac{\sin{\bar{\Omega}}}{\sin{i}}\frac{\partial\mathcal{R}}{\partial\omega}\\\
\end{split}$ (21)

Then, after some algebra, the equation can be written as

$\begin{split}&\frac{d\bar{\Omega}}{dt}\Bigg{|}_{l=2}=\frac{\mathcal{G}M_{\star}^{2}R_{p}^{5}}{a^{6}}\sum_{m=0}^{2}\frac{(2-m)!}{(2+m)!}(2-\delta_{0,m})\Bigg{\\{}\\\
&\Big{(}\frac{1}{C\dot{\theta}\tan{i}}-\frac{\cos{\bar{\Omega}}}{C\dot{\theta}\tan{\varepsilon}}+\frac{1}{\beta
na^{2}\sqrt{1-e^{2}}\sin{i}}\Big{)}\times\\\
&\frac{1}{2}\sum_{p=0}^{2}\frac{\partial\leavevmode\nobreak\
F_{2mp}(i)^{2}}{\partial\leavevmode\nobreak\
i}\sum_{q=-7}^{7}G_{2pq}(e)^{2}\Re\big{(}k_{2}(\sigma_{2mpq})\big{)}\\\
&-\frac{\sin{\bar{\Omega}}\cot{i}}{C\dot{\theta}\tan{\varepsilon}}m\sum_{p=0}^{2}F_{2mp}(i)^{2}\sum_{q=-7}^{7}G_{2pq}(e)^{2}\Im\big{(}k_{2}(\sigma_{2mpq})\big{)}\\\
&+\frac{1}{C\dot{\theta}\tan{\varepsilon}}\frac{\sin{\bar{\Omega}}}{\sin{i}}\sum_{p=0}^{l}(2-2p)F_{2mp}(i)^{2}\sum_{q=-7}^{7}G_{2pq}(e)^{2}\Im\big{(}k_{2}(\sigma_{2mpq})\big{)}\Bigg{\\}}\\\
\end{split}$ (22)

The equation of the argument of the periastron derivative is defined as

$\begin{split}\frac{d\omega}{dt}\Bigg{|}_{l=2}&=-\frac{\beta}{C\dot{\theta}\sin{i}}\frac{\partial\mathcal{R}}{\partial
i}+\frac{\sqrt{1-e^{2}}}{na^{2}e}\frac{\partial\mathcal{R}}{\partial
e}-\frac{\cos{i}}{na^{2}\sqrt{1-e^{2}}\sin{i}}\frac{\partial\mathcal{R}}{\partial
i}\\\ \end{split}$ (23)

Then, the equation can be written as

$\begin{split}&\frac{d\omega}{dt}\Bigg{|}_{l=2}=\frac{\mathcal{G}M_{\star}^{2}R_{p}^{5}}{a^{6}}\sum_{m=0}^{2}\frac{(2-m)!}{(2+m)!}\frac{(2-\delta_{0,m})}{2}\Bigg{[}\\\
&-\Bigg{(}\frac{1}{C\dot{\theta}\sin{i}}+\frac{1}{na^{2}\sqrt{1-e^{2}}\tan{i}}\frac{1}{\beta}\Bigg{)}\sum_{p=0}^{2}\frac{\partial\leavevmode\nobreak\
F_{2mp}(i)^{2}}{\partial\leavevmode\nobreak\
i}\sum_{q=-7}^{7}G_{2pq}(e)\Re\big{(}k_{2}(\sigma_{2mpq})\big{)}\\\
&+\frac{\sqrt{1-e^{2}}}{na^{2}e}\frac{1}{\beta}\sum_{p=0}^{2}F_{2mp}(i)\sum_{q=-7}^{7}\frac{\partial\leavevmode\nobreak\
G_{2pq}(e)^{2}}{\partial\leavevmode\nobreak\
e}\Re\big{(}k_{2}(\sigma_{2mpq})\big{)}\Bigg{]}\\\ \end{split}$ (24)

Finally, the equation of the inclination of the spin axis is defined as

$\begin{split}\frac{d\varepsilon}{dt}&=-\frac{\beta}{C\dot{\theta}}\Big{(}\cos{\bar{\Omega}}\cot{i}\frac{\partial\mathcal{R}}{\partial\bar{\Omega}}+\sin{\bar{\Omega}}\frac{\partial\mathcal{R}}{\partial
i}-\frac{\cos{\bar{\Omega}}}{\sin{i}}\frac{\partial\mathcal{R}}{\partial\omega}\Big{)}\\\
\end{split}$ (25)

Then, after some algebra, we find

$\begin{split}\frac{d\varepsilon}{dt}&=\frac{\mathcal{G}M_{\star}^{2}R^{5}}{a^{6}C\dot{\theta}}\sum_{m=0}^{2}\frac{(2-m)!}{(2+m)!}(2-\delta_{0,m})\Bigg{(}\\\
&m\cos{\bar{\Omega}}\cot{i}\sum_{p=0}^{2}F_{2mp}(i)^{2}\sum_{q=-2}^{2}G_{2pq}(e)^{2}\Im\big{(}k_{2}(\sigma_{2mpq})\big{)}\\\
&-\frac{1}{2}\sin{\bar{\Omega}}\sum_{p=0}^{2}\frac{\partial
F_{2mp}(i)^{2}}{\partial
i}\sum_{q=-\infty}^{\infty}G_{2pq}(e)^{2}\Re\big{(}k_{2}(\sigma_{2mpq})\big{)}\\\
&-\frac{\cos{\bar{\Omega}}}{\sin{i}}\sum_{p=0}^{2}(2-2p)F_{2mp}(i)^{2}\sum_{q=-2}^{2}G_{2pq}(e)^{2}\Im\big{(}k_{2}(\sigma_{2mpq})\big{)}\Bigg{)}\\\
\end{split}$ (26)

##  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 3 and 4. The eccentricity functions are elliptic
expansions that can be computed with the Hansen function
$\mathbf{X}_{l-2p+q}^{-(l-1),(l-2p)}$ (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$.

Table 3: Inclination polynomials $F_{lmp}(i)$ for $l=2$ (Kaula (1964))  l | m | p | $F_{lmp}(i)$
---|---|---|---
2 | 0 | 0 | $-\frac{3}{8}\sin{i}^{2}$
2 | 0 | 1 | $\frac{3}{4}\sin{i}^{2}-\frac{1}{2}$
2 | 0 | 2 | $-\frac{3}{8}\sin{i}^{2}$
2 | 1 | 0 | $\frac{3}{4}\sin{i}(1+\cos{i})$
2 | 1 | 1 | $-\frac{3}{2}\sin{i}\cos{i}$
2 | 1 | 2 | $\frac{3}{4}\sin{i}(\cos{i}-1)$
2 | 2 | 0 | $\frac{3}{4}(1+\cos{i})^{2}$
2 | 2 | 1 | $\frac{3}{2}\sin{i}^{2}$
2 | 2 | 2 | $\frac{3}{4}(1-\cos{i})^{2}$
Table 4: Eccentricity polynomials $G_{lpq}(e)$ from Cayley (1861), up to order $7$ in eccentricity. l | p | q | $G_{lpq}(e)$
---|---|---|---
2 | 0 | -7 | $\frac{15625.0}{129024.0}e^{7}$
2 | 0 | -6 | $\frac{4.0}{45.0}e^{6}$
2 | 0 | -5 | $\frac{81.0}{1280.0}e^{5}+\frac{81.0}{2048.0}e^{7}$
2 | 0 | -4 | $\frac{1.0}{24.0}e^{4}+\frac{7.0}{240.0}e^{6}$
2 | 0 | -3 | $\frac{1.0}{48.0}e^{3}+\frac{11.0}{768.0}e^{5}+\frac{313.0}{30720.0}e^{7}$
2 | 0 | -2 | $0$
2 | 0 | -1 | $-\frac{1}{2}e+\frac{1}{16}e^{3}-\frac{5}{384}e^{5}-\frac{143}{18432}e^{7}$
2 | 0 | 0 | $1-\frac{5}{2}e^{2}+\frac{13}{16}e^{4}-\frac{35}{288}e^{6}$
2 | 0 | 1 | $\frac{7}{2}e-\frac{123}{16}e^{3}+\frac{489}{128}e^{5}-\frac{1763}{2048}e^{7}$
2 | 0 | 2 | $\frac{17}{2}e^{2}-\frac{115}{16}e^{4}+\frac{601}{48}e^{6}$
2 | 0 | 3 | $\frac{845.0}{48.0}e^{3}-\frac{32525.0}{768.0}e^{5}+\frac{208225.0}{6144.0}e^{7}$
2 | 0 | 4 | $\frac{533.0}{16.0}e^{4}-\frac{13827.0}{160.0}e^{6}$
2 | 0 | 5 | $\frac{228347.0}{3840.0}e^{5}-\frac{3071075.0}{18432.0}e^{7}$
2 | 0 | 6 | $\frac{73369.0}{720.0}e^{6}$
2 | 0 | 7 | $\frac{12144273.0}{71680.0}e^{7}$
2 | 1 | -7 | $G_{217}(e)$
2 | 1 | -6 | $G_{216}(e)$
2 | 1 | -5 | $G_{215}(e)$
2 | 1 | -4 | $G_{214}(e)$
2 | 1 | -3 | $G_{213}(e)$
2 | 1 | -2 | $G_{212}(e)$
2 | 1 | -1 | $G_{211}(e)$
2 | 1 | 0 | $(1-e^{2})^{-3/2}\simeq 1+\frac{3}{2}e^{2}+\frac{15}{8}e^{4}+\frac{35}{16}e^{6}+\mathcal{O}(e^{9})$
2 | 1 | 1 | $\frac{3.0}{2.0}e+\frac{27.0}{16.0}e^{3}+\frac{261.0}{128.0}e^{5}+\frac{14309.0}{6144.0}e^{7}$
2 | 1 | 2 | $\frac{9.0}{4.0}e^{2}+\frac{7.0}{4.0}e^{4}+\frac{141.0}{64.0}e^{6}$
2 | 1 | 3 | $\frac{53.0}{16.0}e^{3}+\frac{393.0}{256.0}e^{5}+\frac{24753.0}{10240.0}e^{7}$
2 | 1 | 4 | $\frac{77.0}{16.0}e^{4}+\frac{129.0}{160.0}e^{6}$
2 | 1 | 5 | $\frac{1773.0}{256.0}e^{5}-\frac{4987.0}{6144.0}e^{7}$
2 | 1 | 6 | $\frac{3167.0}{320.0}e^{6}$
2 | 1 | 7 | $\frac{432091.0}{30720.0}e^{7}$
2 | 2 | -7 | $G_{207}(e)$
2 | 2 | -6 | $G_{206}(e)$
2 | 2 | -5 | $G_{205}(e)$
2 | 2 | -4 | $G_{204}(e)$
2 | 2 | -3 | $G_{203}(e)$
2 | 2 | -2 | $G_{202}(e)$
2 | 2 | -1 | $G_{201}(e)$
2 | 2 | 0 | $G_{200}(e)$
2 | 2 | 1 | $G_{20-1}(e)$
2 | 2 | 2 | $G_{20-2}(e)$
2 | 2 | 3 | $G_{20-3}(e)$
2 | 2 | 4 | $G_{20-4}(e)$
2 | 2 | 5 | $G_{20-5}(e)$
2 | 2 | 6 | $G_{20-6}(e)$
2 | 2 | 7 | $G_{20-7}(e)$