∎

11institutetext: C. Ragazzo (ORCID 0000-0002-4277-4173)22institutetext:
Instituto de Matemática e Estatística, Universidade de São Paulo, 05508-090
São Paulo, SP, Brazil
22email<EMAIL_ADDRESS>
Partially supported by FAPESP grant 2016/25053-8. 33institutetext: L. S. Ruiz
(ORCID 0000-0002-5705-5278) 44institutetext: Instituto de Matemática e
Computação, Universidade Federal de Itajubá, 37500-903 Itajubá, MG, Brazil
44email<EMAIL_ADDRESS>
# Tidal Evolution and Spin-Orbit Dynamics: The Critical Role of Rheology.

Clodoaldo Ragazzo  Lucas Ruiz dos Santos

###### Abstract

This study analyzes secular dynamics using averaged equations that detail
tidal effects on the motion of two extended bodies in Keplerian orbits. It
introduces formulas for energy dissipation within each body of a binary
system. The equations, particularly in contexts like the Sun-Mercury system,
can be delineated into a fast-slow system. A significant contribution of this
work is the demonstration of the crucial role complex rheological models play
in the capture by spin-orbit resonances. This is particularly evident in the
notable enlargement of the basin of attraction for Mercury’s current state
when transitioning from a single characteristic time rheology to a dual
characteristic time model, under the constraint that both models comply with
the same estimate of the complex Love number at orbital frequency. The study
also underscores the importance of Mercury’s elastic rigidity on secular
timescales.

##  1 Introduction

Newton formulated the law of gravitation and concluded that the motion of the
centers of mass of spherical bodies is equivalent to that of point masses. The
solution to the resulting two-body problem, as obtained by Newton, forms the
backbone of all subsequent developments in celestial mechanics. Notably, while
planets and major satellites are almost spherical, even slight deformations
caused by spin and tidal forces can significantly influence their rotation and
orbits.

Building on the foundational works of Newton, Laplace, Thomson, Darwin, and
others, a significant advancement in understanding tidal effects on celestial
motion was made by Kaula Kaula (1964), who decomposed tidal forces
harmonically in space and time for two bodies in a Keplerian orbit, using Love
numbers to evaluate changes in orbital elements. Recent updates to Kaula’s
theory are available in Boué and Efroimsky (2019), with additional insights in
Efroimsky (2012). Over the last 70 years, extensive research has focused on
these tidal effects. Notable contributions include Ferraz-Mello’s Ferraz-Mello
(2013) model, which builds on Darwin’s theory with non-spherical hydrostatic
states, detailed further in Folonier et al. (2018), Ferraz-Mello et al.
(2020), and summarized in Ferraz-Mello (2019), Ferraz-Mello (2021). Various
studies, including those by Goldreich (1966), Singer (1968), Alexander (1973),
Mignard (1979), Hut (1981), Makarov and Efroimsky (2013), Correia et al.
(2014), Ferraz-Mello (2015), and Boué et al. (2016), have explored deformation
equations averaged over orbital motion, particularly in low and high-viscosity
scenarios. The averaged equations used in this work are identical to those
presented in Correia and Valente (2022).

In this paper we study the secular-planar dynamics of two extended bodies. We
make the following assumptions:

  * $1)$

The two bodies are deformable, nearly spherical at all times;

  * $2)$

The spins (or rotation vectors) of the deformable bodies remain perpendicular
to the orbital plane.

  * 3)

The bodies are: radially stratified, each body layer is homogeneous and has a
linear visco-elastic rheology, see e.g. Sabadini et al. (2016). Fluid layers,
if present, must be sufficiently coupled to the adjacent layers such that the
rotation of each layer remains close to the average rotation of the body.

In this paper, we crucially use that, from the perspective of gravitation, the
rheology of a body with a finite number of homogeneous layers is equivalent to
that of a homogeneous body with a sufficiently more complex rheology Gevorgyan
et al. (2023).

The aim of this study is to analyze the secular dynamics arising from the
averaged equations. The primary novelty of this work lies in demonstrating the
importance of using complex rheological models in the capture by spin-orbit
resonances. More specifically, it highlights the significant enlargement of
the basin of attraction of Mercury’s current state when transitioning from a
rheology with one characteristic time to one with two characteristic times.
Both rheologies are constrained to satisfy the same estimate of the complex
Love number at orbital frequency.

In the next section, we present the main results of the paper. The final
section provides several mathematical details involved in deriving these
results.

##  2 Main results.

The foundational equations for the orbit and rotation of an extended body are
well-established in the literature. Various equations detailing the
deformation of extended bodies exist. In a companion paper Ragazzo and Ruiz
(2024), we averaged the equations provided in Ragazzo and Ruiz (2017) with
respect to the orbital motion, excluding the term accounting for the inertia
of deformations as discussed in Correia et al. (2018). These equations are
applicable to any rheological model. We obtained essentially the same averaged
equations as those presented in Correia and Valente (2022). The only
difference is the inclusion of a centrifugal deformation term, which is not
relevant in the planar case. It is important to emphasize that the averaged
equations in Correia and Valente (2022) are more general as they do not
necessitate the spins of the bodies to be perpendicular to the orbital plane.

Let $m_{\alpha}$ and $m$ represent the masses of two celestial bodies, which
could be a planet and a star, or a planet and a satellite, etc. We name the
bodies such that the “$\alpha$ body” is the largest $m_{\alpha}\geq m$. The
quantities of the large body will always be labeled with an index $\alpha$ and
those of the small body will have no label.

We assume that both bodies are almost spherical, deformable, and the
deformations are volume preserving. In this situation the mean moments of
inertia, denoted as ${\,\rm I}_{\circ}$ and ${\rm I}_{\circ\alpha}$, remain
constant in time, a result attributed to Darwin Rochester and Smylie (1974).

The motion of the two-body system is determined by three sets of equations:
the equations for the relative positions of the centers of mass (orbit), the
equations for the rotation of each body about their centers of mass (spin),
and the equations for the deformation of each body. In Section 3.1, we present
these fundamental equations.

###  2.1 Rheology and Love numbers

The rheology of a body determines the Love number $k_{2}(\sigma)$, where
$\sigma$ is the angular frequency of the tidal raising force. For a stratified
body with a finite number of homegeneous layers, the Love number can be
written as Sabadini et al. (2016) (see also Gevorgyan et al. (2023) and
Gevorgyan (2021)):

$\begin{split}k_{2}(\sigma)&=k_{\infty}+(k_{\circ}-k_{\infty})\left(\frac{h_{1}}{1+i\tau_{1}\sigma}+\cdots+\frac{h_{n+1}}{1+i\tau_{n+1}\sigma}\right).\
\ \end{split}$ (2.1)

In these equations:

  * $\bullet$

$k_{\infty}:=\lim_{\sigma\to\infty}k_{2}(\sigma)$ is the asymptotic value of
the Love number at high frequencies.

  * $\bullet$

$k_{\circ}:=k_{2}(0)$ is the Love number at frequency zero, also called
secular Love number 111 The softest possible body is one composed of a perfect
fluid, which is held together solely by self-gravity. In this case,
$k_{\circ}=k_{f}$, where $k_{f}$ is the fluid Love number. Let $R$ be the body
volumetric radius and $R_{I}$ be the radius of inertia, defined as the radius
of a homogeneous sphere of mass $m$ and moment of inertia ${\,\rm I}_{\circ}$.
Assuming that the density is non-decreasing towards the center, the following
approximation (valid for $0.2<\frac{{\,\rm I}_{\circ}}{mR^{2}}\leq 0.4$) holds
(Ragazzo, 2020, Eqs. 1.2 and 1.8, Theorems 4.1 and 4.2):
$k_{f}\approx\frac{3}{2}\left(\frac{R_{I}}{R}\right)^{5}\quad\text{where:}\quad
R_{I}^{2}:=\frac{5}{2}\frac{{\,\rm I}_{\circ}}{m}\,.$ (2.2) The maximum value
of $R_{I}/R$ is one, achieved in a homogeneous body, for which
$k_{f}=\frac{3}{2}$. The approximation in equation (2.2) was also proposed in
Consorzi et al. (2023)..

  * $\bullet$

$\tau_{i}$ is the $i^{th}$ characteristic time of the rheology.

  * $\bullet$

$h_{i}>0$ is the relative amplitude of the $i^{th}$ mode of the rheology,
which is associated to the characteristic time $\tau_{i}$. The relative
amplitudes add to one:

$h_{1}+h_{2}+\ldots h_{n+1}=1.$

As discussed in Gevorgyan et al. (2023), the Love number in equation (2.1)
corresponds to that of a homogeneous body with a generalized Voigt rheology,
the spring-dashpot representation of which is given in Figure 1.

Figure 1: The generalized Voigt model. The constant $\gamma$ (blue) represents
self-gravity rigidity.

All the simple rheological models used in the literature, such as: Maxwell,
Kelvin-Voigt, Standard Anelastic Solid (SAS), Burgers, Bland’s generalized
Voigt (Figure 2 LEFT), and Bland’s generalized Maxwell (Figure 2 RIGHT); are
particular cases of the generalized Voigt rheology depicted in Figure 1. The
Andrade model can be approximated with arbitrary precision by a generalized
Voigt model, as demonstrated in Gevorgyan et al. (2020).

Figure 2:  LEFT: Generalized Voigt model in Bland (1960) Figure 1.8. RIGHT:
Generalized Maxwell model Bland (1960) Figure 1.9.

###  2.2 Passive deformation.

The moment of inertia matrix $\mathbf{I}$ of a body in an inertial reference
frame can be written as:

$\mathbf{I}={\,\rm I}_{\circ}\big{(}{\bf 1}-\mathbf{b}\big{)}.$ (2.3)

where ${\bf 1}$ is the identity and $\mathbf{b}$ is a symmetric and traceless
matrix. We denote matrices and vectors in bold face. The matrix $\mathbf{b}$
is termed the deformation matrix. The deformation matrix of the $\alpha$-body
is denoted as $\mathbf{b}_{\alpha}$.

Since the bodies are nearly spherical and tidal deformations are minor, the
variations in orbital elements occur gradually. As a result, the tidal forces
can be approximated by those of point masses undergoing Keplerian motion.
Within this framework, tidal forces can be harmonically decomposed both
temporally, using Hansen coefficients222 If $a$ is the semi-major axis, $e$ is
the eccentricity, and $M$ is the mean anomaly, then
$\left(\frac{r}{a}\right)^{n}{\rm
e}^{imf}=\sum_{k=-\infty}^{\infty}X^{n,m}_{k}(e){\rm e}^{ikM}\,,$ (2.4) where:
$n$, $m$, and $k$ are integers, and $X^{n,m}_{k}(e)$ is the Hansen
coefficient., and spatially, using spherical harmonics.

In Section 3.2, we employ these decompositions to compute the “passive
deformation matrix” in terms of Love numbers, orbital eccentricity, mean
motion ($n=\dot{M}$), and the spin angular velocities ($\omega$ and
$\omega_{\alpha}$) of each body. The term “passive” refers to the fact that
this deformation is computed while neglecting its influence on altering the
orbital elements and the spin of each body.

###  2.3 Energy and angular momentum.

The energy function of a system of two rigid spherical bodies is

$E_{\circ}:=-\frac{Gmm_{\alpha}}{2a}+{\,\rm
I}_{\circ}\frac{\omega^{2}}{2}+{\,\rm
I}_{\circ\alpha}\frac{\omega_{\alpha}^{2}}{2}=-a_{1}n^{2/3}+{\,\rm
I}_{\circ}\frac{\omega^{2}}{2}+{\,\rm
I}_{\circ\alpha}\frac{\omega_{\alpha}^{2}}{2}\,,$ (2.5)

where: $G$ is the gravitational constant, $a$ is the semi-major axis, and

$a_{1}:=\frac{mm_{\alpha}G^{2/3}}{(m+m_{\alpha})^{1/3}}\,.$ (2.6)

The angular momentum of a spherical body is $\ell_{s}={\,\rm I}_{\circ}\omega$
and the orbital angular momentum is

$\ell=\frac{mm_{\alpha}}{m+m_{\alpha}}a^{2}n\sqrt{1-e^{2}}=a_{1}n^{-1/3}\sqrt{1-e^{2}}\,.$
(2.7)

The total angular momentum of a system of two spherical bodies is

$\ell_{\scriptscriptstyle T}:=\ell+\ell_{s}+\ell_{s\alpha}.$ (2.8)

We will assume that $\ell>0$ (or $n>0$) and $\ell_{\scriptscriptstyle T}>0$
for all time.

The energy and angular momentum of a system of two slightly deformable bodies,
which are given in Section 3.1, are approximately given by the expressions of
their spherical approximations.

The total angular momentum is constant, $\dot{\ell}_{\scriptscriptstyle T}=0$,
and within the Lagrangian formalism with dissipation function Ragazzo and Ruiz
(2017) the time derivative of the energy is given by the sum of the
dissipation functions of each body $\cal D$ and $\cal D_{\alpha}$.

In Section 3.3 we use the passive deformations to estimate the average
dissipation of energy in each body. For the small body the result is:

$\begin{split}\langle\mathcal{D}\rangle&=-\frac{3}{8}\left(\frac{m_{\alpha}}{m+m_{\alpha}}\right)^{2}\frac{n^{4}R^{5}}{G}\sum_{k=-\infty}^{\infty}\Bigg{\\{}\frac{kn}{3}\Big{(}X^{-3,0}_{k}(e)\Big{)}^{2}{\rm
Im}\,k_{2}(kn)\\\ &\ +(kn-2\omega)\Big{(}X^{-3,2}_{k}(e)\Big{)}^{2}{\rm
Im}\,k_{2}(kn-2\omega)\Bigg{\\}}.\end{split}$ (2.9)

A similar expression holds for the $\alpha$-body.

Since Im $k_{2}(-\sigma)=-$Im $k_{2}(\sigma)$ and Im $k_{2}(\sigma)<0$ for
$\sigma>0$, then $\langle\mathcal{D}\rangle\geq 0$. If $e=0$, then
$X^{-3,0}_{k}(0)=0$ for $k\in{\mathbb{Z}}$ and and $X^{-3,2}_{k}(0)=0$ for
$k\neq 2$. So $\langle\mathcal{D}\rangle=0$ if and only if $e=0$ and
$n=\omega$.

The energy dissipated due to passive deformations must originate from the
motion of the spherical bodies that induce these passive deformations.
Consequently, we can infer

$\begin{split}\dot{E}_{\circ}&=-\frac{2}{3}a_{1}n^{-1/3}\dot{n}+{\,\rm
I}_{\circ}\omega\dot{\omega}+{\,\rm
I}_{\circ\alpha}\omega_{\alpha}\dot{\omega}_{\alpha}\\\
&=-2\langle\mathcal{D}\rangle-2\langle\mathcal{D_{\alpha}}\rangle.\end{split}$
(2.10)

###  2.4 The average torque and the secular equation for the orbital elements

For the small body, the average torque due to passive deformations is:

$\begin{split}\langle{\rm
T}\rangle&=-\frac{3}{2}\left(\frac{m_{\alpha}}{m+m_{\alpha}}\right)^{2}\frac{n^{4}R^{5}}{G}\bigg{\\{}\sum\limits_{k=-\infty}^{\infty}\Big{(}X^{-3,2}_{k}\Big{)}^{2}{\rm
Im}\,k_{2}(kn-2\omega)\bigg{\\}}.\end{split}$ (2.11)

A similar expression holds for the $\alpha$-body.

The secular equation for the orbital elements can be readily derived from
equations (2.9), (2.10), (2.11), and their counterparts for the $\alpha$-body:

$\begin{split}\dot{\omega}&=\frac{1}{{\,\rm I}_{\circ}}\langle{\rm
T}\rangle,\\\ \dot{\omega}_{\alpha}&=\frac{1}{{\,\rm
I}_{\circ\alpha}}\langle{\rm T}_{\alpha}\rangle,\\\
\dot{n}&=\frac{3n^{1/3}}{2a_{1}}\Big{(}2\langle\mathcal{D}\rangle+2\langle\mathcal{D_{\alpha}}\rangle+\omega\langle{\rm
T}\rangle+\omega_{\alpha}\langle{\rm T}_{\alpha}\rangle\Big{)}.\end{split}$
(2.12)

The conservation of total angular momentum, combined with equations (2.7) and
(2.8), implies that the eccentricity, which is featured on the right-hand side
of equations (2.12), can be expressed in terms of the state variables
$\omega,\omega_{\alpha},n$.

###  2.5 Time scales and a simplification when $m\ll m_{\alpha}$.

The despin rate is dependent on the imaginary parts of the Love numbers. The
largest bodies in the solar system are fluid (e.g., the Sun, Jupiter, etc.)
and have an imaginary part of the Love number that is significantly smaller
than that of bodies with solid components. Therefore, it is reasonable to
assume in the subsequent equations that either ${\rm Im}\,k_{2}$ and ${\rm
Im}\,k_{2\alpha}$ are comparable or ${\rm Im}\,k_{2}\gg{\rm Im}\,k_{2\alpha}$.

The ratio between the despin rates, assuming no spin-orbit resonances, is
given by

$\frac{\dot{\omega}}{\dot{\omega}_{\alpha}}\approx\frac{{\,\rm
I}_{\circ\alpha}}{{\,\rm
I}_{\circ}}\frac{m_{\alpha}^{2}}{m^{2}}\frac{R^{5}}{R_{\alpha}^{5}}\frac{{\rm
Im}\,k_{2}}{{\rm
Im}\,k_{2\alpha}}\approx\frac{\rho_{\alpha}}{\rho}\frac{m_{\alpha}^{2}}{m^{2}}\frac{{\rm
Im}\,k_{2}}{{\rm Im}\,k_{2\alpha}}\,,$ (2.13)

where $\rho$ represents the density of the body. If $m_{\alpha}\gg m$, then
$\dot{\omega}\gg\dot{\omega}_{\alpha}$, indicating that the despin rate of the
larger body is significantly slower than that of the smaller body. In such a
case, as a first approximation, we may assume $\omega_{\alpha}=$constant and
equations (2.12) become:

$\begin{split}\dot{\omega}&=\frac{1}{{\,\rm I}_{\circ}}\langle{\rm
T}\rangle,\\\
\dot{n}&=\frac{3n^{1/3}}{2a_{1}}\Big{(}2\langle\mathcal{D}\rangle+2\langle\mathcal{D_{\alpha}}\rangle+\omega\langle{\rm
T}\rangle\Big{)}.\end{split}$ (2.14)

For example, for the Sun-Mercury system, using Im$k_{2\alpha}\approx 3.5\times
10^{-8}$ Ogilvie (2014) and Im$k_{2}=0.0051$ Margot et al. (2018), we obtain
$\frac{\dot{\omega}}{\dot{\omega}_{\alpha}}\approx 2.8\times 10^{17}$.

###  2.6 A simplification when the $\alpha$ body is fluid (almost inviscid).

The ratio of the energy dissipation rates in each body is given by

$\frac{\langle{\cal D}\rangle}{\langle{\cal
D}_{\alpha}\rangle}\approx\frac{m_{\alpha}^{2}}{m^{2}}\frac{R^{5}}{R_{\alpha}^{5}}\frac{{\rm
Im}\,k_{2}}{{\rm
Im}\,k_{2\alpha}}\approx\frac{\rho_{\alpha}^{2}}{\rho^{2}}\frac{R_{\alpha}}{R}\frac{{\rm
Im}\,k_{2}}{{\rm Im}\,k_{2\alpha}}.$ (2.15)

Equation (2.15), along with the equation for $\dot{n}$ in (2.14), implies that
depending on the imaginary parts of the Love numbers, the larger body may play
a more significant role than the smaller body in altering the orbital elements
333For the Earth-Moon system, using the Love number
$k_{2\alpha}=0.2817-0.02324\,i$ for the Earth (Ragazzo and Ruiz, 2017, Table
3, semi-diurnal frequency), and ${\rm Im}\,k_{2}=5.13\times 10^{-4}$ for the
Moon (orbital frequency) Fienga et al. (2019), we obtain  $\frac{\langle{\cal
D}\rangle}{\langle{\cal
D}_{\alpha}\rangle}\approx\frac{m_{\alpha}^{2}R^{5}\,{\rm
Im}\,k_{2}}{m^{2}R_{\alpha}^{5}\,{\rm Im}\,k_{2\alpha}}=0.022.$ (2.16) This
indicates that the Earth’s influence on altering the orbital parameters is
greater than that of the Moon. .

However, when the $\alpha$ body consists of a low viscosity fluid and the
smaller body includes a solid part, then $\frac{{\rm Im}\,k_{2}}{{\rm
Im}\,k_{2\alpha}}\gg 1$, and $\frac{\langle{\cal D}\rangle}{\langle{\cal
D}_{\alpha}\rangle}$ may also be significantly greater than one. This scenario
applies to the Sun-Mercury system, where $\frac{\langle{\cal
D}\rangle}{\langle{\cal D}_{\alpha}\rangle}\approx 2.8\times 10^{6}$. In such
cases, $\langle{\cal D}_{\alpha}\rangle$ can be neglected in comparison to
$\langle{\cal D}\rangle$ in the equation for $\dot{n}$ in (2.14).

Therefore, for the Sun-Mercury system, we can employ the following
approximation to equation (2.12):

$\begin{split}\dot{\omega}&=\frac{1}{{\,\rm I}_{\circ}}\langle{\rm
T}\rangle,\\\
\dot{n}&=\frac{3n^{1/3}}{2a_{1}}\Big{(}2\langle\mathcal{D}\rangle+\omega\langle{\rm
T}\rangle\Big{)}\,,\end{split}$ (2.17)

where the state variables are $(\omega,n)$. This type of equation commonly
appears in the literature, as seen in references such as Correia et al.
(2014), Gomes et al. (2019), and Correia and Valente (2022).

The angular momentum, which remains constant along the motion, becomes
simpler:

$\ell_{\scriptscriptstyle T}=a_{1}n^{-1/3}\sqrt{1-e^{2}}+{\,\rm
I}_{\circ}\omega.$ (2.18)

Equations (2.17) are equivalent to those obtained in Correia and Valente
(2022) in the planar case.

###  2.7 The dynamics when $\frac{{\,\rm I}_{\circ}}{ma^{2}}\ll 1$.

In the remainder of this section, we consider the simplest case where equation
(2.17) is applicable, with $m\ll m_{\alpha}$ and ${\rm Im}\,k_{2\alpha}\ll{\rm
Im}\,k_{2}$.

The ratio

$\begin{split}\frac{\dot{n}}{\dot{\omega}}&=\frac{{\,\rm
I}_{\circ}n^{4/3}}{a_{1}}\frac{3\omega}{2n}\left\\{1+{\scriptstyle\frac{\frac{1}{2\omega}\Bigg{\\{}\sum\limits_{k=-\infty}^{\infty}\frac{kn}{3}\Big{(}X^{-3,0}_{k}(e)\Big{)}^{2}{\rm
Im}\,k_{2}(kn)+\ldots\Bigg{\\}}}{\sum\limits_{k=-\infty}^{\infty}\Big{(}X^{-3,2}_{k}\Big{)}^{2}{\rm
Im}\,k_{2}(kn-2\omega)}}\right\\}\\\ \end{split}$ (2.19)

includes the factor

$\frac{{\,\rm
I}_{\circ}n^{4/3}}{a_{1}}=\frac{1}{m}\left(1+\frac{m}{m_{\alpha}}\right)\frac{{\,\rm
I}_{\circ}}{a^{2}}\approx\frac{{\,\rm I}_{\circ}}{ma^{2}}\,.$ (2.20)

We assume

$\frac{{\,\rm I}_{\circ}n^{4/3}}{a_{1}}\approx\frac{{\,\rm
I}_{\circ}}{ma^{2}}\ll 1\,.$ (2.21)

For the Sun-Mercury system,

$\frac{{\,\rm I}_{\circ}n^{4/3}}{a_{1}}=6.2\times 10^{-10}\,.$ (2.22)

If $\left|\frac{\omega}{n}\right|$ is not far from one and $(n,\omega)$ is not
on

${\cal
C}:=\bigg{\\{}(n,\omega):\sum\limits_{k=-\infty}^{\infty}\Big{(}X^{-3,2}_{k}(e)\Big{)}^{2}{\rm
Im}\,k_{2}(kn-2\omega)=0\bigg{\\}}\quad\text{\bf(torque-free curve)}\,,$
(2.23)

then $\frac{\dot{\omega}}{\dot{n}}\gg 1$. As a result, the spin-orbit dynamics
is governed by a slow-fast system of equations, where $n$ is the slow variable
and $\omega$ is the fast variable. The typical dynamics is described in the
caption of Figure 3.

Figure 3: Typical dynamics when $\frac{ma^{2}}{{\,\rm I}_{\circ}}\gg 1$. We
use the nondimensional mean motion $\nu=\frac{\ell^{3}_{\scriptscriptstyle
T}}{a^{3}_{1}}\,n$ and the nondimensional spin angular velocity ${\rm
w}=\frac{\ell^{3}_{\scriptscriptstyle T}}{a^{3}_{1}}\,\omega$ in the plot, as
seen in equation (2.26) and the previous paragraphs. The curve in black
represents the torque-free curve $\cal C$: $\dot{\rm w}<0$ above $\cal C$,
$\dot{\rm w}>0$ below $\cal C$, and $\dot{\rm w}=0$ on $\cal C$. An initial
condition above $\cal C$ converges rapidly to $\cal C$ with a decreasing spin
velocity, while an initial condition below $\cal C$ converges rapidly to $\cal
C$ with an increasing spin velocity. Upon reaching $\cal C$, the solution
gradually progresses along $\cal C$. The decrease in energy implies that $\nu$
increases along $\cal C$ (with $e$ decreasing). The motion along $\cal C$ is
sustainable up to a fold point. At this point, the decrease in energy forces
the solution to depart from $\cal C$, resulting in a spin jump. The solution
is then drawn towards another point on $\cal C$, representing a lower-order
spin-orbit resonance. Eventually, the solution approaches the synchronous
solution $(\nu,{\rm w})\approx(1,1)$. For a general discussion on the flow
close to $\mathcal{C}$, see Ragazzo and Ruiz (2024).

Due to inequality (2.21), equation (2.18) can be simplified to

$\frac{\ell_{\scriptscriptstyle T}}{n{\,\rm I}_{\circ}}=\frac{a_{1}}{{\,\rm
I}_{\circ}n^{4/3}}\sqrt{1-e^{2}}+\frac{\omega}{n}\approx\frac{a_{1}}{{\,\rm
I}_{\circ}n^{4/3}}\sqrt{1-e^{2}}\,,$ (2.24)

where we have used that $\left|\frac{\omega}{n}\right|$ is not far from one
and $e$ is not close to one. This leads to the approximation

$e\approx\sqrt{1-\left(\frac{\ell^{3}_{\scriptscriptstyle
T}}{a^{3}_{1}}n\right)^{2/3}}.$ (2.25)

The constant $\frac{\ell^{3}_{\scriptscriptstyle T}}{a^{3}_{1}}$ has the unit
of time (for the Sun-Mercury system $\frac{\ell^{3}_{\scriptscriptstyle
T}}{a^{3}_{1}}=13.1224$ days444 $2\pi\frac{\ell^{3}_{\scriptscriptstyle
T}}{a^{3}_{1}}=82.45$ days is approximately the period of the longest possible
synchronous, $\omega=n$, circular orbit of Mercury about the Sun.). We will
use this constant to define the nondimensional mean motion $\nu$ and spin
angular velocity ${\rm w}$, which are used in several figures:

$\nu:=\frac{\ell^{3}_{\scriptscriptstyle T}}{a^{3}_{1}}\,n\qquad{\rm
w}:=\frac{\ell^{3}_{\scriptscriptstyle T}}{a^{3}_{1}}\,\omega.$ (2.26)

###  2.8 Characteristic times of terrestrial planets.

The imaginary part of the Love number of Mercury remains largely unknown, even
at its orbital frequency Baland et al. (2017), Steinbrügge et al. (2018). The
most well-studied terrestrial bodies in are the Earth and the Moon. It is
reasonable to hypothesize that the rheology of Mercury could be similar, at
least qualitatively, to those of the Earth and Moon.

Thirteen characteristic times and amplitudes of a five-layer interior model
for the Moon Matsuyama et al. (2016) were computed in Gevorgyan et al. (2023).
This research demonstrated that the complex Love number of the model can be
effectively represented within the frequency range of interest by just three
modes. These modes have nondimensional amplitudes $h$ and characteristic times
$\tau$ (in days) (Gevorgyan et al. (2023) Table 2 555 The correspondence
between the notation $(s_{j},r_{j})$ in Gevorgyan et al. (2023) and that used
here is: $s_{j}=-\tau_{j}^{-1}$, $r_{j}=h_{j}(k_{\circ}-k_{\infty})$. The
values of $h_{1},h_{2},h_{3}$ in (2.27) are normalized such that
$h_{1}+h_{2}+h_{3}=1$):

$\begin{split}(h_{1},\tau_{1})&=\big{(}0.487,12.2\,{\rm d}\big{)}\,,\\\
(h_{2},\tau_{2})&=\big{(}0.512,18.3\,{\rm d}\big{)}\,,\\\
(h_{3},\tau_{3})&=\big{(}0.001,8.3\times 10^{6}\,{\rm d}\big{)}.\\\
\end{split}$ (2.27)

A five-layer interior model for the Earth is presented in Table 2-1 in
Sabadini et al. (2016). Nine relaxation times for this model, ranging from
$2.6$ to $2\times 10^{5}$ centuries, are provided in (Sabadini et al., 2016,
p.64), but without corresponding amplitudes.

These examples demonstrate that terrestrial bodies, such as Mercury:

  * $\bullet$

May possess more than one relevant characteristic time,

  * $\bullet$

The relevant characteristic times can differ by orders of magnitude.

To investigate the importance of Mercury’s rheology in the spin-orbit
dynamics, we compare two different scenarios. In the first, Mercury is assumed
to have a rheology with exactly one characteristic time $\tau$. In the second
scenario, Mercury is assumed to have a rheology with exactly two
characteristic times $\tau_{1}$ and $\tau_{2}$, such that
$\tau_{1}\ll\tau_{2}$ and $\tau_{2}n_{mer}\gg 1$, where $n_{mer}$ is the
orbital mean motion.

###  2.9 Mercury’s Love Number.

For Mercury, at the orbital frequency

$n_{mer}=\frac{2\pi}{87.969\text{ days}}\,,$ (2.28)

current estimates of the Love number fall within the range $0.53<{\rm
Re}\,k_{2}<0.63$, with the preferred value being $k_{2}=0.569$ Genova et al.
(2019), Goossens et al. (2022), and $0<-{\rm Im}\,k_{2}<0.025$ (Baland et al.,
2017, p. 152), (Steinbrügge et al., 2018, p.2767 and 2769). We will assume
${\rm Re}\,k_{2}=0.57$ and $-{\rm Im}\,k_{2}=0.011$ as reference values
(quality factor $Q\approx-\frac{{\rm Re}\,k_{2}}{{\rm Im}\,k_{2}}\approx 52$).

Mercury is not in hydrostatic equilibrium: the expected hydrostatic flattening
is two orders of magnitude lower than the observed flattening. This
discrepancy in hydrostaticity may be attributed to Mercury’s despinning
history Matsuyama and Nimmo (2009). According to Ragazzo et al. (2022)
(Section 4.2), this suggests the presence of secular-elastic rigidity,
represented by $\alpha_{0}>0$ in Figure 1, in conjunction with self-gravity,
such that $k_{\circ}<k_{f}$, where $k_{f}$ denotes the fluid Love number.
Equation (2.2) implies that $k_{f}=1.04$.

The choice of $k_{\circ}$ within the range $[0.57,1.04]$ is critical. One of
the main results we present below – the significant enlargement of the basin
of attraction of Mercury’s current state when transitioning from a rheology
with one characteristic time to one with two characteristic times – does not
hold if $k_{\circ}=k_{f}$. This highlights the importance of Mercury’s secular
elastic rigidity in determining its current state. In our subsequent analysis,
we will use $k_{\circ}=0.8$ and $k_{\circ}=0.7$.

###  2.10  Mercury with One Relaxation Time.

The generalized Voigt rheology (or standard solid rheology) with a single
characteristic time is represented by a spring-dashpot model shown in Figure 1
with $n=0$. The Love number is described by equation (2.1) with $n=0$:

$k_{2}(\sigma)=k_{\infty}+(k_{\circ}-k_{\infty})\,.\frac{1}{1+i\tau\sigma}.\\\
$ (2.29)

In Figure 4, we show the torque-free curve for the choice $k_{\circ}=0.8$.
Imposing $k_{2}(n_{mer})=0.57-0.011\,i$ we obtain $k_{\infty}=0.5695$ and
$\tau n_{mer}=20.91$. Both the current state of Mercury and the stable
synchronous equilibrium are also depicted. The current state of Mercury, as
illustrated in Figure 4, suggests that Mercury must have been captured
directly into the 3-2 spin-orbit resonance, without prior transitions through
higher-order resonances.

In Figure 5, we show the torque-free curve for the choice $k_{\circ}=0.7$.
Imposing $k_{2}(n_{mer})=0.57-0.011\,i$ we obtain $k_{\infty}=0.5691$ and
$\tau n_{mer}=11.89$. In this scenario, the current state of Mercury suggests
that it could have been captured in a higher-order resonance before settling
into the 3-2 spin-orbit resonance. The potential initial states of Mercury
that lead to its current state, as indicated by the region hatched with
vertical lines in Figure 5, encompass a much larger set than the corresponding
region in Figure 4.

Figure 4: Torque-free curve (blue) for Mercury with $k_{\circ}=0.8$ and one
relaxation time: $\tau n_{merc}=20.91$. Three spin-orbit resonance lines
(dashed) are shown. Mercury is currently on the $3n=2\omega$ resonance line.
The red point represents the synchronous equilibrium. The region hatched with
vertical black lines represents potential initial positions of Mercury that
could have led to its current state. The variables in the horizontal and
vertical axes are: $\nu=\frac{\ell^{3}_{\scriptscriptstyle T}}{a^{3}_{1}}\,n$,
${\rm w}=\frac{\ell^{3}_{\scriptscriptstyle T}}{a^{3}_{1}}\,\omega$, where
$\frac{\ell^{3}_{\scriptscriptstyle T}}{a^{3}_{1}}=\frac{0.9373}{n_{mer}}$.
Figure 5: Torque-free curve (blue) for Mercury with $k_{\circ}=0.7$ and one
relaxation time $\tau n_{mer}=11.82$. The remainder of the notation is the
same as that in Figure 4.

###  2.11 Mercury with Two Relaxation Times.

We assume a generalized Voigt rheology, equation (2.1) and Figure 1 with
$n=1$.

We choose $k_{\circ}=0.8$ and fix the complex Love number at the orbital
frequency as $k_{2}(n_{\text{mer}})=0.57-0.011\,i$. Equation (2.1) then
implies

$0.57-0.011\,i=k_{\infty}+(0.8-k_{\infty})\left(\frac{h_{1}}{1+i\tau_{1}n_{\text{mer}}}+\frac{h_{2}}{1+i\tau_{2}n_{\text{mer}}}\right)$
(2.30)

with $h_{1}>0$, $h_{2}>0$, and $h_{1}+h_{2}=1$. This results in a set of two
scalar equations for the four unknowns: $k_{\infty}$, $0<h_{1}<1$,
$\tau_{1}n_{\text{mer}}$, and $\tau_{2}n_{\text{mer}}$.

In Section 3.4, we demonstrate that equation (2.30) has solutions with
arbitrarily large values of $\frac{\tau_{2}}{\tau_{1}}$. According to the
models discussed in Section 2.8, the ratio between the largest and smallest
relaxation times of both the Moon and Earth is of the order of $10^{5}$.
Therefore, it would be logical to choose solutions for equation (2.30) where
$\frac{\tau_{2}}{\tau_{1}}$ is at least of order $10^{5}$. However, this
choice presents numerical challenges in plotting the torque-free curve.
Consequently, some plots presented below have a smaller ratio of
$\frac{\tau_{2}}{\tau_{1}}$. As argued in Section 3.4, the conclusions drawn
would remain valid for any larger value of $\frac{\tau_{2}}{\tau_{1}}$.

In Figure 6, we present the torque-free curve for the values:
$\tau_{1}n_{\text{mer}}=0.48,\tau_{2}n_{\text{mer}}=1396,h_{1}=0.11,k_{\infty}=0.5474$,
which implies $k_{2}(n_{\text{mer}})=0.57-0.11\,i$. In this scenario, the
current state of Mercury suggests that it could have initially been captured
in a higher-order resonance before transitioning to the 3-2 spin-orbit
resonance. The potential initial states of Mercury that could lead to its
current state, as indicated by the region hatched with vertical lines in
Figure 6, encompass a much larger set than the corresponding region in Figure
4, which was derived using a rheology with a single relaxation time under the
same constraints $k_{2}(n_{\text{mer}})=0.57-0.11\,i$ and $k_{\circ}=0.8$.

Figure 6: Torque-free curve (blue) for Mercury with:
$\tau_{1}n_{mer}=0.48,\tau_{2}n_{mer}=1396,h_{1}=0.11,h_{2}=0.89,k_{\infty}=0.5474$.
The remainder of the notation is the same as that in Figure 4.

In Figure 7, we illustrate the effect on the torque-free curve of increasing
the ratio $\frac{\tau_{2}}{\tau_{1}}$ by a factor of 10. In both cases, the
current state of Mercury suggests that it might have initially been captured
in a higher-order resonance before transitioning to the 3-2 spin-orbit
resonance. The conclusions regarding the set of initial conditions that lead
to Mercury’s current state, as established in Figure 6, are equally applicable
to Figure 7.

Figure 7: Torque-free curve (blue) for Mercury. LEFT:
$\tau_{1}n_{\text{mer}}=0.015,\tau_{2}n_{\text{mer}}=94,h_{1}=0.7126,k_{\infty}=0$;
and RIGHT:
$\tau_{1}n_{\text{mer}}=0.019,\tau_{2}n_{\text{mer}}=940,h_{1}=0.7128,k_{\infty}=0$.

The parameters used to generate Figure 6 do not accurately reflect the current
state of Mercury, where $\omega=\frac{3}{2}n_{\text{mer}}$. According to
equation (2.17) and with the parameters employed in Figure 6, the current rate
of variation of Mercury’s spin is $\frac{\dot{\omega}}{\omega}=\frac{1}{{\,\rm
I}_{\circ}\omega}\langle{\rm T}\rangle=-2.5\times 10^{-7}\,\text{yr}^{-1}$. At
this rate, over 10 million years, Mercury’s spin would decrease to $10\%$ of
its current value. The same conclusion applies to the parameters used to
generate Figure 7 and many other set of parameters that we tested. This may
indicate that trapping into the $3:2$ spin-orbit resonance happened by means
of a spin decrease and not the opposite.

The strong decay rate of the spin seems paradoxical in light of Mercury’s
current position in Figure 6, which appears to be on the torque-free curve.
This paradox is elucidated in Figure 8 left, where we can observe that the
gradient of the torque function at a point on the torque-free curve is very
large. Consequently, the torque may still be nonnegligible even at points
quite close to the torque-free curve. In Figure 8 right, we demonstrate that
for $n=n_{\text{mer}}$, the value of $\frac{\omega}{n_{\text{mer}}}$ where the
torque is zero is slightly below $\frac{3}{2}$, so Mercury is positioned above
the torque-free curve.

Figure 8: Two graphs of $\dot{\omega}=\frac{\langle T\rangle}{{\,\rm
I}_{\circ}}$ divided by the constant
$(k_{\circ}-k_{\infty})\frac{3n_{\text{mer}}^{4}R^{5}}{2G{\,\rm
I}_{\circ}}=3.33\times 10^{-19}\,\text{s}^{-2}$ are presented as a function of
$\frac{\omega}{n_{\text{mer}}}$. The parameters are those used in Figure 6.
The graph on the left displays the two roots of $\dot{\omega}$, one of which
appears to be at $\frac{\omega}{n_{\text{mer}}}=1.5$. The graph on the right
demonstrates that the root close to $\frac{\omega}{n_{\text{mer}}}=1.5$ is
actually slightly less than this value.

###  2.12 Conclusion

In this paper, we have obtained, by means of energy arguments, equations for
the secular evolution of spin and orbit due to tidal effects. These equations
were originally obtained by other authors, notably in Correia and Valente
(2022). In the case of Mercury, which served as a model in our study on the
dynamical effects of increasingly complex rheologies, the equations can be
further simplified.

The main conclusion, which may be applicable to any binary system and not just
Mercury-Sun, is that employing more complex and realistic rheologies is
crucial in understanding spin-orbit evolution. This conclusion aligns with the
findings in Makarov (2012) and Noyelles et al. (2014). According to the
abstract of Noyelles et al. (2014): “As opposed to the commonly used constant
time lag (CTL) and constant phase lag (CPL) models, the physics-based tidal
model changes dramatically the statistics of the possible final spin states.”
Here, the term “physics-based tidal model” refers to a model derived from a
rheology with several relaxation times.

The past capture of Mercury into its current spin-orbit state has been the
subject of several previous studies. In our paper, we have completely
neglected several important effects that may have been significant, as pointed
out in Noyelles et al. (2014): perturbations by other bodies in the Solar
System that induce variations in Mercury’s eccentricity, the potential
relevance of core-mantle-boundary friction at fluid-solid interfaces, the
possibility of past impacts with other bodies, existence of a permanent
deformation, and changes in the rheology over long time scales. Within this
broader context, two aspects of Mercury’s rheology, which were crucial in our
model, may still remain significant: 1) The elastic rigidity of Mercury at
secular time scales, which must exist; otherwise, Mercury would be in
hydrostatic equilibrium. 2) The presence of two or more relaxation times in
Mercury’s rheology. Both these factors dramatically alter the probability of
Mercury being captured in its current state. However, our quantitative results
regarding the spin-orbit history of Mercury almost certainly do not align with
what actually happened.

##  3 Mathematical Complements

In this section, we present several mathematical facts that explain and
clarify many of the statements made in the previous section.

###  3.1 The fundamental equations

Let $\kappa:=\\{\mathbf{e}_{1},\mathbf{e}_{2},\mathbf{e}_{3}\\}$ be an
inertial orthonormal frame with origin at the centre of mass of the system.
The orbit is contained in the plane $\\{\mathbf{e}_{1},\mathbf{e}_{2}\\}$.

Let $\mathbf{x}$ be the vector from the center of mass of the small body to
the body $\alpha$. The equation for the orbital motion is

$\begin{split}\ddot{\mathbf{x}}&=G(m_{\alpha}+m)\bigg{\\{}-\frac{\mathbf{x}}{|\mathbf{x}|^{3}}-\frac{15}{2}\frac{1}{|x|^{7}}\left(\left(\frac{{\,\rm
I}_{\circ\alpha}}{m_{\alpha}}\mathbf{b}_{\alpha}+\frac{{\,\rm
I}_{\circ}}{m}\mathbf{b}\right)\mathbf{x}\cdot\mathbf{x}\right)\mathbf{x}\\\
&+3\frac{1}{|\mathbf{x}|^{5}}\left(\frac{{\,\rm
I}_{\circ\alpha}}{m_{\alpha}}\mathbf{b}_{\alpha}+\frac{{\,\rm
I}_{\circ}}{m}\mathbf{b}\right)\mathbf{x}\bigg{\\}}\,,\end{split}$ (3.31)

where the deformation matrix of the small body $\mathbf{b}$, defined in
equation (2.3), is given by

$\mathbf{b}=\left(\begin{array}[]{ccc}b_{11}&b_{12}&0\\\ b_{12}&b_{22}&0\\\
0&0&b_{33}\\\ \end{array}\right)\,,\quad\text{with}\quad
b_{33}=-b_{11}-b_{22}\,.$ (3.32)

For a rigid body, there exists a frame, namely the body frame, in which the
body remains at rest. Specifically, the angular momentum of the body with
respect to the body frame is null. In the case of a deformable body, a frame
still exists where the body’s angular momentum is null: this is known as the
Tisserand frame. Let $\mathrm{K}:=\\{\mathbf{e}_{{\scriptscriptstyle
T1}},\mathbf{e}_{{\scriptscriptstyle T2}},\mathbf{e}_{{\scriptscriptstyle
T3}}\\}$ denote an orthonormal Tisserand frame for the small body. The
orientation of $\mathrm{K}$ with respect to the inertial frame
$\kappa=\\{\mathbf{e}_{1},\mathbf{e}_{2},\mathbf{e}_{3}\\}$ is described by
$\mathbf{R}:\mathrm{K}\to\kappa$. We assume that the angular velocity of the
small body, $\boldsymbol{\omega}$, is perpendicular to the orbital plane,
which implies :

$\mathbf{R}(\phi)=\left(\begin{array}[]{ccc}\cos\phi&-\sin\phi&0\\\
\sin\phi&\cos\phi&0\\\ 0&0&1\end{array}\right)\,,$ (3.33)

with $\boldsymbol{\omega}=\omega\mathbf{e}_{3}=\dot{\phi}\mathbf{e}_{3}$.

The angular momentum of the small body is denoted by
$\boldsymbol{\ell}_{s}=\ell_{s}\mathbf{e}_{3}$, with the index $s$
representing spin, and is given by:

$\ell_{s}=\omega{\,\rm I}_{\circ}(1-b_{33})\,.$ (3.34)

In the context of the quadrupolar approximation, Euler’s equation for the
variation of $\ell_{s}$ is:

$\dot{\ell}_{s}=-\frac{3G{\,\rm
I}_{\circ}m_{\alpha}}{\|\mathbf{x}\|^{5}}\Big{\\{}x_{1}x_{2}(b_{22}-b_{11})+b_{12}\left(x_{1}^{2}-x_{2}^{2}\right)\Big{\\}}.$
(3.35)

A similar equation holds for the large body.

It is imperative to note that, unlike the case of a rigid body, the
deformation matrix is not constant in the Tisserand frame $\mathrm{K}$.

To complete the set of equations (3.31) and (3.35), we require additional
equations for the deformation matrices. These equations were derived within
the Lagrangian formalism and utilizing what was termed the “Association
Principle”, as detailed in Ragazzo and Ruiz (2015), Ragazzo and Ruiz (2017).

In the Tisserand frame $\mathrm{K}$ of the small body, the deformation matrix
and the position vector are denoted by capital letters as follows:

$\mathbf{B}=\mathbf{R}(\phi)\mathbf{b}\mathbf{R}^{-1}(\phi)\qquad\mathbf{X}=\mathbf{R}^{-1}(\phi)\mathbf{x}\,.$
(3.36)

The equations for the deformation variables of the smaller body for the
generalised Voigt model in Figure 1 are (Ragazzo et al., 2022, Eq. (3.15))

$\begin{split}&(\gamma+\alpha_{0})\mathbf{B}+\mathbf{\Lambda}=\mathbf{F}\\\
&\frac{1}{\alpha}\dot{\mathbf{\Lambda}}+\frac{1}{\eta}\mathbf{\Lambda}=\dot{\tilde{\mathbf{B}}}\\\
&\eta_{j}\dot{\mathbf{B}}_{j}+\alpha_{j}\mathbf{B}_{j}=\mathbf{\Lambda},\quad
j=1,\ldots,n\\\
&\mathbf{B}=\tilde{\mathbf{B}}+\mathbf{B}_{1}+\mathbf{B}_{2}+\ldots+\mathbf{B}_{n}\,,\end{split}$
(3.37)

where:

  * $\bullet$

$\gamma$, with dimensions of $1/$time2, is a parameter representing the self-
gravity rigidity of the body; a larger $\gamma$ indicates a stronger
gravitational force holding the body together. In Figure 1 $\gamma$ would be a
spring in parallel to the spring-dashpot system representing the rheology.

  * $\bullet$

$\alpha_{0}$, $\alpha$, $\alpha_{1},\ldots$ also with dimensions of $1/$time2,
are elastic rigidity coefficients; $\alpha$ increases with the stiffness.

  * $\bullet$

$\eta$, $\eta_{1},\ldots$, dimensions of $1/$time, are viscosity parameters; a
body with a larger $\eta$ is harder to deform at a given rate compared to a
body with a smaller $\eta$.

  * $\bullet$

$\tilde{\mathbf{B}}$, $\mathbf{B}_{1},\ldots$, are nondimensional traceless
matrices that represent internal variables of the rheology;

  * $\bullet$

$\mathbf{\Lambda}$, with dimensions 1/time2, is a traceless force matrix that
represents the stress upon the rheology part ($\mathbf{F}-\mathbf{\Lambda}$ is
the stress supported by self-gravity).

  * $\bullet$

$\mathbf{F}$, with dimensions 1/time2, is the force matrix in the Tisserand
frame $\mathrm{K}$:

$\begin{array}[]{l l l}\mathbf{F}&:=\mathbf{C}+\mathbf{S}&\text{Deformation
force}\\\ &&\\\ \mathbf{C}&:=\frac{\omega^{2}}{3}\left(\begin{matrix}1&0&\ \
0\\\ 0&1&\ \ 0\\\ 0&0&-2\end{matrix}\right)&\text{Centrifugal force}\\\ &&\\\
\mathbf{S}&:=\frac{3Gm_{\alpha}}{|\mathbf{X}|^{5}}\left(\mathbf{X}\mathbf{X}^{T}-\frac{|\mathbf{X}|^{2}}{3}{\bf
1}\right)&\text{Tidal force}\end{array}$ (3.38)

where $\mathbf{X}$ is undestood as a column vector, $\mathbf{X}^{T}$ is the
transpose of $\mathbf{X}$ (a row vector vector), and
$\mathbf{X}\mathbf{X}^{T}$ is the usual product of a $n\times 1$ by $1\times
n$ matrix, which gives an $n\times n$ matrix with entries
$\big{(}\mathbf{X}\mathbf{X}^{T}\big{)}_{ij}=X_{i}X_{j}$. The norm of a matrix
is defined as $\|\mathbf{B}\|^{2}=\frac{1}{2}{\rm\\!\ Tr\\!\
}\big{(}\mathbf{B}^{T}\mathbf{B}\big{)}$.

Conservation of angular momentum. The orbital angular momentum is given by

$\ell\mathbf{e}_{3}=\boldsymbol{\ell}=\frac{m_{\alpha}m}{m_{\alpha}+m}\mathbf{x}\times\dot{\mathbf{x}}\,.$
(3.39)

From equations (3.31) and (3.35) we obtain

$\dot{\boldsymbol{\ell}}=3\frac{m_{\alpha}m}{|\mathbf{x}|^{5}}\mathbf{x}\times\left(\frac{{\rm
I}_{\circ\alpha}}{m_{\alpha}}\mathbf{B}_{\alpha}+\frac{{\,\rm
I}_{\circ}}{m}\mathbf{B}\right)\mathbf{x}=-\dot{\boldsymbol{\ell}}_{s\alpha}-\dot{\boldsymbol{\ell}}_{s}\,,$
(3.40)

that shows the conservation of total angular momentum

$\boldsymbol{\ell}_{{\scriptscriptstyle
T}}:=\boldsymbol{\ell}+\boldsymbol{\ell}_{s\alpha}+\boldsymbol{\ell}_{s}\,.$
(3.41)

Energy. The total energy of the system is

$E=E_{kcm}+E_{krot}+E_{def}+E_{gr}$ (3.42)

where:

$E_{kcm}=\frac{m_{\alpha}m}{m_{\alpha}+m}\frac{|\dot{\mathbf{x}}|^{2}}{2}$
(3.43)

is the kinectic energy of the orbital motion;

$E_{krot}=\frac{\boldsymbol{\omega}_{\alpha}\cdot\mathbf{I}_{\alpha}\boldsymbol{\omega}_{\alpha}}{2}+\frac{\boldsymbol{\omega}\cdot\mathbf{I}\boldsymbol{\omega}}{2}$
(3.44)

is the kinectic energy of rotation;

$E_{def}=\frac{{\,\rm
I}_{\circ}}{2}\left((\gamma+\alpha_{0})\|\mathbf{B}\|^{2}+\frac{\|\mathbf{\Lambda}\|^{2}}{\alpha}+\sum_{j=1}^{n}\alpha_{j}\|\mathbf{B}_{j}\|^{2}\right)+\
\text{Large body term}$ (3.45)

is the elastic energy of deformation of the generalised Voigt model in Figure
1; and

$\begin{split}E_{gr}=&-\frac{Gm_{\alpha}m}{|\mathbf{x}|}-\frac{3G}{2}\frac{1}{|\mathbf{x}|^{5}}\bigl{\\{}m{\rm
I}_{\circ\alpha}\big{(}\mathbf{x}\cdot\mathbf{B}_{\alpha}\mathbf{x}\big{)}+m_{\alpha}{\,\rm
I}_{\circ}(\mathbf{x}\cdot\mathbf{B}\mathbf{x})\bigr{\\}}\end{split}$ (3.46)

is the potential energy due to gravitational interactions.

Dissipation of energy. For the generalised Voigt model in Figure 1, the
dissipation function of the small body is:

$\mathcal{D}=\frac{{\,\rm
I}_{\circ}}{2}\left(\frac{\|\mathbf{\Lambda}\|^{2}}{\eta}+\sum_{j=1}^{n}\eta_{j}\|\dot{\mathbf{B}}_{j}\|^{2}\right)$
(3.47)

A similar expression $\mathcal{D}_{\alpha}$ holds for the large body. The
Lagrangian formalism with dissipation function Ragazzo and Ruiz (2015),
Ragazzo and Ruiz (2017) implies that the total energy decreases along the
motion:

$\dot{E}=-(2\mathcal{D}+2\mathcal{D}_{\alpha})\leq 0\,.$ (3.48)

Moreover, $2\mathcal{D}_{\alpha}$ and $2\mathcal{D}$ are the powers dissipated
within the large and small bodies, respectively.

The secular equations (2.12) are approximations to the fundamental equations
presented in this Section.

###  3.2 Computation of the Passive Deformation Matrix.

Passive deformation refers to the tides caused by an orbiting point mass on an
extended body, while the influence of the gravitational field, resulting from
the deformation on the orbit, is neglected. The equations that describe these
passive deformations are obtained setting $\mathbf{B}=\mathbf{B}_{\alpha}=0$
in equations (3.31) and (3.35) while preserving equation (3.37):

$\begin{split}&\ddot{\mathbf{x}}=-G(m_{\alpha}+m)\frac{\mathbf{x}}{|\mathbf{x}|^{3}}\,,\quad\dot{\omega}=0\,,\quad\dot{\omega}_{\alpha}=0\\\
&(\gamma+\alpha_{0})\mathbf{B}+\mathbf{\Lambda}=\mathbf{F}\\\
&\frac{1}{\alpha}\dot{\mathbf{\Lambda}}+\frac{1}{\eta}\mathbf{\Lambda}=\dot{\tilde{\mathbf{B}}}\\\
&\eta_{j}\dot{\mathbf{B}}_{j}+\alpha_{j}\mathbf{B}_{j}=\mathbf{\Lambda},\quad
j=1,\ldots,n\\\
&\mathbf{B}_{T}=\tilde{\mathbf{B}}+\mathbf{B}_{1}+\mathbf{B}_{2}+\ldots+\mathbf{B}_{n}\,.\\\
&\text{Plus deformation equations for the large body}\end{split}$ (3.49)

In this scenario $\omega$ and $\omega_{\alpha}$ remain constant and
$\mathbf{x}$ follows a Keplerian ellipse.

In order to write the $\mathbf{F}$ in a convenient way, we define the matrices
(Ragazzo and Ruiz, 2017, Eq. (41))

$\mathbf{Y}_{0}:=\frac{1}{\sqrt{3}}\left(\begin{matrix}1&0&\ \ 0\\\ 0&1&\ \
0\\\
0&0&-2\end{matrix}\right)\,,\quad\mathbf{Y}_{1}:=\frac{1}{\sqrt{2}}\left(\begin{matrix}0&0&1\\\
0&0&-i\\\
1&-i&0\end{matrix}\right)\,,\quad\mathbf{Y}_{2}:=\frac{1}{\sqrt{2}}\left(\begin{matrix}1&-i&0\\\
-i&-1&0\\\ 0&0&0\end{matrix}\right)\,,$ (3.50)

$\mathbf{Y}_{-1}=\overline{\mathbf{Y}}_{1}$, and
$\mathbf{Y}_{-2}=\overline{\mathbf{Y}}_{2}$, where the overline represents
complex conjugation.

Any symmetric matrix can be writen as a linear combination of the six matrices
$\\{{\bf
1},\mathbf{Y}_{-2},\mathbf{Y}_{-1},\mathbf{Y}_{0},\mathbf{Y}_{1},\mathbf{Y}_{2}\\}$.
This basis is orthonormal with respect to the Hermitean inner product
$\mathbf{A}\cdot\mathbf{B}=\frac{1}{2}{\rm\\!\ Tr\\!\
}\big{(}\,\overline{\mathbf{A}}^{\rm
T}\mathbf{B}\big{)}=\frac{1}{2}\sum_{ij}\overline{A}_{ij}B_{ij}$.

Traceless matrices that are invariant under rotation with respect to the
$\mathbf{e}_{3}$-axis can be written as linear combinations of
$\\{\mathbf{Y}_{-2},\mathbf{Y}_{0},\mathbf{Y}_{2}\\}$. These matrices have a
simple transformation rule with respect to rotations about the axis
$\mathbf{e}_{3}$, namely

$\mathbf{R}(\theta)\mathbf{Y}_{j}\mathbf{R}^{-1}(\theta)={\rm
e}^{i\,j\,\theta}\,\mathbf{Y}_{j}\,,\quad j=-2,0,2\,.$ (3.51)

An elliptic orbit in the body frame is given by

$\begin{split}\mathbf{X}&=\mathbf{R}^{-1}(\phi)\mathbf{x}=r\mathbf{R}(f+\varpi-\phi)\mathbf{e}_{1}\\\
&=r(\cos(f+\varpi-\phi)\,\mathbf{e}_{1}+\sin(f+\varpi-\phi)\,\mathbf{e}_{2})\,,\end{split}$
(3.52)

where $f$ is the mean anomaly, $\varpi$ is the argument of the periapsis.

The associated tidal-force matrix, equation (3.38), can be written as

$\begin{split}\mathbf{S}&=\frac{3Gm_{\alpha}}{r^{3}}\mathbf{R}_{3}(f+\varpi-\phi)\Big{\\{}\mathbf{e}_{1}\mathbf{e}_{1}^{T}-\frac{1}{3}{\bf
1}\Big{\\}}\mathbf{R}^{-1}_{3}(f+\varpi-\phi)\\\
&=\frac{3Gm_{\alpha}}{2r^{3}}\bigg{\\{}{\rm
e}^{-2i(f+\varpi-\phi)}\frac{\mathbf{Y}_{-2}}{\sqrt{2}}+\frac{\mathbf{Y}_{0}}{\sqrt{3}}+{\rm
e}^{2i(f+\varpi-\phi)}\frac{\mathbf{Y}_{2}}{\sqrt{2}}\bigg{\\}}.\end{split}$
(3.53)

In equation (3.53), the variables $r$, $f$, and $\phi=\omega t$ are dependent
on $t$.

Equations (3.53) and (2.4) gives a harmonic decomposition of the tidal force
as follows:

$\mathbf{S}=\frac{3Gm_{\alpha}}{2a^{3}}\sum_{l=-2}^{2}\sum_{k=-\infty}^{\infty}{\rm
e}^{i\\{t(kn-l\omega)+l\,\varpi\\}}\mathbf{Y}_{l}U_{kl},$ (3.54)

where $U_{k,-1}=U_{k,1}=0$ and

$U_{k,-2}=\frac{X^{-3,-2}_{k}}{\sqrt{2}}\,,\ \
U_{k0}=\frac{X^{-3,0}_{k}}{\sqrt{3}}\,,\ \
U_{k2}=\frac{X^{-3,2}_{k}}{\sqrt{2}}\,.$ (3.55)

The symmetry property $X^{n,-m}_{-k}=X^{n,m}_{k}$ implies $U_{kj}=U_{-k,-j}$.

The centrifugal force in equation (3.38) can be represented as

$\mathbf{C}=\frac{\omega^{2}}{\sqrt{3}}\mathbf{Y}_{0}\,.$ (3.56)

To obtain the almost periodic solution of the deformation equation solving for
each Fourier mode separately suffices.

To simplify the notation we consider a simple harmonic force term of the form

$\mathbf{F}(t)=\mathbf{F}^{\prime}{\rm e}^{i\,\sigma\,t}$

where $\mathbf{F}^{\prime}$ is a complex amplitude matrix, and
$\sigma\in{\mathbb{R}}$ is the constant forcing frequency.

If we do the substitutions

$\mathbf{B}\to\mathbf{B}^{\prime}\,{\rm
e}^{i\,\sigma\,t}\,,\quad\mathbf{B}_{j}\to\mathbf{B}_{j}^{\prime}\,{\rm
e}^{i\,\sigma\,t}\,,\quad\mathbf{\Lambda}\to\mathbf{\Lambda}^{\prime}\,{\rm
e}^{i\,\sigma\,t}\,,$

where $\mathbf{B}^{\prime}$, $\mathbf{B}_{j}^{\prime}$ and
$\mathbf{\Lambda_{j}}^{\prime}$ are understood as constant complex-amplitude
matrices, into equation (3.37), then we obtain after some simplifications

$\left\\{\gamma+\alpha_{0}+\left(\frac{1}{\alpha}+\frac{1}{\eta
i\sigma}+\sum_{j=1}^{n}\frac{1}{\alpha_{j}+i\sigma\eta_{j}}\right)^{-1}\right\\}\mathbf{B}^{\prime}=\mathbf{F}^{\prime}\,.$
(3.57)

The term $\frac{1}{\alpha}+\frac{1}{\eta
i\sigma}+\sum_{j=1}^{n}\frac{1}{\alpha_{j}+i\sigma\eta_{j}}$ is the complex
rigidity $J(\sigma)$ of the generalised Voigt model when $\alpha_{0}=0$
(Ragazzo et al., 2022, equation 4.23).

The complex Love number at frequency $\sigma$ can be written as (Correia et
al., 2018, Section 4):

$k_{2}(\sigma)=\frac{3{\,\rm
I}_{\circ}G}{R^{5}}\frac{1}{\gamma+\alpha_{0}+J^{-1}(\sigma)}\,.$ (3.58)

Note: $k_{2}(-\sigma)=\overline{k}_{2}(\sigma)$. Therefore, equation (3.57)
can be written as

$\mathbf{B}^{\prime}=\frac{R^{5}}{3{\,\rm
I}_{\circ}G}k_{2}(\sigma)\mathbf{F}^{\prime}.$ (3.59)

Applying equation (3.59) to each Fourier coefficient of
$\mathbf{F}=\mathbf{S}+\mathbf{C}$, where $\mathbf{S}$ and $\mathbf{C}$ are
defined in equations (3.54) and (3.56) respectively, we obtain the passive-
deformation matrix

$\begin{split}\mathbf{B}(t)&=k_{\circ}\zeta_{c}\frac{\mathbf{Y}_{0}}{3\sqrt{3}}+\zeta_{\scriptscriptstyle
T}\sum_{l=-2}^{2}\sum_{k=-\infty}^{\infty}{\rm
e}^{i\\{t(kn-l\omega)+l\,\varpi\\}}k_{2}(kn-l\omega)\mathbf{Y}_{l}U_{kl},\end{split}$
(3.60)

where $k_{\circ}=k_{2}(0)$ is the secular Love number and

$\zeta_{c}:=\frac{R^{5}\omega^{2}}{G{\,\rm
I}_{\circ}}\quad\text{and}\quad\zeta_{\scriptscriptstyle
T}:=\frac{m_{\alpha}R^{5}}{2{\,\rm I}_{\circ}a^{3}}.$ (3.61)

###  3.3 The Average Dissipation of Energy and Average Torque.

The energy dissipated in the smaller body is $-2\mathcal{D}$, where
$\mathcal{D}$ is the dissipation function as given in equation (3.47). The
average dissipation function over the orbital motion is approximately given by

$\langle\mathcal{D}\rangle:=\lim_{t\to\infty}\int\limits_{0}^{T}\mathcal{D}dt=\frac{{\,\rm
I}_{\circ}}{2}\lim_{t\to\infty}\int\limits_{0}^{T}\left(\frac{\|\mathbf{\Lambda}\|^{2}}{\eta}+\sum_{j=1}^{n}\eta_{j}\|\dot{\mathbf{B}}_{j}\|^{2}\right)dt\,,$
(3.62)

where $(\mathbf{B},\mathbf{\Lambda},\ldots)$ are the components of the
solution to equation (3.37) under the passive deformation hypothesis.

From equation (3.37) we obtain:

$\begin{split}&(\gamma+\alpha_{0})\bcancel{\langle\dot{\mathbf{B}}\cdot\mathbf{B}\rangle}+\langle\dot{\mathbf{B}}\cdot\mathbf{\Lambda}\rangle=\bcancel{\langle\dot{\mathbf{B}}\cdot\mathbf{C}\rangle}+\langle\dot{\mathbf{B}}\cdot\mathbf{S}\rangle\\\
&\bcancel{\frac{1}{\alpha}\langle\mathbf{\Lambda}\cdot\dot{\mathbf{\Lambda}}\rangle}+\frac{1}{\eta}\|\mathbf{\Lambda}\|^{2}=\langle\mathbf{\Lambda}\cdot\dot{\tilde{\mathbf{B}}}\rangle\\\
&\eta_{j}\|\dot{\mathbf{B}}_{j}\|^{2}+\alpha_{j}\bcancel{\langle\dot{\mathbf{B}}_{j}\cdot\mathbf{B}_{j}\rangle}=\langle\dot{\mathbf{B}}_{j}\cdot\mathbf{\Lambda}\rangle,\quad
j=1,\ldots,n\\\
&\langle\dot{\mathbf{B}}\cdot\mathbf{\Lambda}\rangle=\langle\dot{\tilde{\mathbf{B}}}\cdot\mathbf{\Lambda}\rangle+\langle\dot{\mathbf{B}}_{1}\cdot\mathbf{\Lambda}\rangle+\langle\dot{\mathbf{B}}_{2}\cdot\mathbf{\Lambda}\rangle+\ldots+\langle\dot{\mathbf{B}}_{n}\cdot\mathbf{\Lambda}\rangle\,,\end{split}$
(3.63)

where we used that the centrifugal force $\mathbf{C}$, equation (3.56), is
constant. These equations imply

$\langle\mathcal{D}\rangle=\frac{{\,\rm
I}_{\circ}}{2}\langle\dot{\mathbf{B}}\cdot\mathbf{S}\rangle\,,$ (3.64)

where $\mathbf{S}$ and $\mathbf{B}$ are given in equations (3.54) and (3.60),
respectively.

In the computation of the averaging
$\langle\dot{\mathbf{B}}\cdot\mathbf{S}\rangle$ the following identities are
useful: $\mathbf{Y}_{l}\cdot\mathbf{Y}_{-m}=\delta_{lm}$, $U_{kl}=U_{-k,-l}$,
and $k_{2}(-\sigma)=\overline{k}_{2}(\sigma)$. After some computations we
obtain equation (2.9).

The total energy, as given in equation (3.42), includes terms that depend on
the deformation matrices. If we set $\mathbf{B}=\mathbf{B}_{\alpha}=0$ in this
equation and assume that the orbit is Keplerian and the spins are constant,
then we obtain the energy function $E_{\circ}$, as presented in equation
(2.5), for a system of two rigid spherical bodies. It is a natural assumption
that the energy dissipated in the two bodies is subtracted from $E_{\circ}$.
This assumption leads to equation (2.10).

The average torque is computed as the average of the right-hand side of
equation (3.35),

$\langle{\rm T}\rangle:=\lim_{t\to\infty}\int\limits_{0}^{T}-\frac{3G{\,\rm
I}_{\circ}m_{\alpha}}{\|\mathbf{x}\|^{5}}\Big{\\{}x_{1}x_{2}(b_{22}-b_{11})+b_{12}\left(x_{1}^{2}-x_{2}^{2}\right)\Big{\\}}dt\,,$
(3.65)

where $\mathbf{B}$ represents the passive deformation matrix and $\mathbf{x}$
determines Keplerian ellipses. Standard computations lead to equation (2.11).

Equations (2.10) and (2.11) imply equations (2.12).

###  3.4 Solutions to Equation (2.30) for $\tau_{2}n_{\text{mer}}\gg 1$ and
the Torque-Free Curve for $\tau_{2}/\tau_{1}\gg 1$

If $\tau_{2}n_{\text{mer}}\gg 1$, then equation (2.30) can be approximately
written as

$\frac{0.57-k_{\infty}}{0.8-k_{\infty}}-\frac{0.011}{0.8-k_{\infty}}\,i=\frac{h_{1}}{1+i\tau_{1}n_{\text{mer}}}-i\frac{1-h_{1}}{\tau_{2}n_{\text{mer}}}.$
(3.66)

The left-hand side defines a curve in the complex plane, parameterized by
$k_{\infty}\in[0,0.57]$. As $\tau_{2}\to\infty$, the right-hand side converges
to $\frac{h_{1}}{1+i\tau_{1}n_{\text{mer}}}$, and for a given $h_{1}\in(0,1)$,
it defines a curve in the complex plane, parameterized by
$\tau_{1}\in(0,\infty)$. This curve is a semi-circle below the real axis, with
a radius of $\frac{h_{1}}{2}$ and centered at
$\left(\frac{h_{1}}{2},0\right)$. In the limit as $\tau_{2}\to\infty$, for a
given $h_{1}$, the solution to equation (3.66) is represented by a pair
$(k_{\infty},\tau_{1})$, corresponding to the intersection of these two
curves, as illustrated in Figure 9.

Figure 9: Graphical solution to equation (3.66) in the limit as
$\tau_{2}\to\infty$ for $h_{1}=0.11$. The dashed-red curve is parameterized by
$k_{\infty}$, while the blue curve is parameterized by $\tau_{1}$. The
intersections of these curves represent the solutions
$(k_{\infty},\tau_{1}n_{\text{mer}})$, specifically: $(0.547658,0.492353)$,
which is marked by a ball, and $(0.563915,1.80765)$.

In the following, we graphically demonstrate how to obtain a solution to
equation (2.30) for $\tau_{2}n_{\text{mer}}\gg 1$, using the solution marked
with a ball in Figure 9.

For the fixed value of $h_{1}=0.11$, we choose $\tau_{1}n_{\text{mer}}$ to be
slightly smaller than $0.492353$, which corresponds to the solution marked in
Figure 9, say $\tau_{1}n_{\text{mer}}=0.48$. The curve parameterized by
$\tau_{1}\in[0,0.48/n_{\text{mer}}]$ (blue) will terminate slightly above the
curve parameterized by $k_{\infty}$ (red). Substituting $h_{1}=0.11$ and
$\tau_{1}n_{\text{mer}}=0.48$ into equation (3.66), we obtain

$\frac{0.57-k_{\infty}}{0.8-k_{\infty}}-\frac{0.011}{0.8-k_{\infty}}\,i=\frac{0.11}{1+i0.48}-i\frac{0.89}{\tau_{2}n_{\text{mer}}}.$
(3.67)

Now, the right-hand side of equation (3.67) defines a vertical segment in the
complex plane starting at $\frac{0.11}{1+i0.48}$ for $\tau_{2}=+\infty$, and
decreasing to $-\infty\,i$ as $\tau_{2}\to 0$. This curve clearly intersects
the curve parameterized by $k_{\infty}$, as depicted in Figure 10. As
$\tau_{1}n_{\text{mer}}$ approaches $0.492353$ from below, the value of
$\tau_{2}n_{\text{mer}}$ increases towards $+\infty$.

An application of the implicit function theorem to equation (3.66) around the
solution
$(k_{\infty},\tau_{1}n_{\text{mer}},(\tau_{2}n_{\text{mer}})^{-1})=(0.547658,0.492353,0)$
yields the following formula:

$k_{\infty}=0.547658-\frac{0.341496}{\tau_{2}n_{mer}}\,,\quad\tau_{1}n_{mer}=0.492353-\frac{17.5779}{\tau_{2}n_{mer}}.$
(3.68)

Figure 10: Graphical solution to equation (3.67) for $\tau_{2}n_{mer}\gg 1$
and $h_{1}=0.11$. The dashed-dotted-black line is parameterized by $\tau_{2}$.

As previously mentioned, plotting the torque-free curve is numerically
challenging when $\tau_{1}\ll\tau_{2}$. In this limit, it is possible to
approximate the torque-free curve near a given resonance line,
$k_{*}n-2\omega=0$, as demonstrated in the following paragraphs.

First, we substitute the variable $\omega$ with $\epsilon$ using the equation
$k_{*}n-2\omega=-2\epsilon$. The analysis will be confined to a neighborhood
of the resonance line, namely

$|\epsilon|\tau_{2}<1\,.$

We write equation (2.23), for the torque-free curve, using the new variable

$\begin{split}0&=-\frac{1}{k_{\circ}-k_{\infty}}\sum\limits_{k=-\infty}^{\infty}\Big{(}X^{-3,2}_{k}(e)\Big{)}^{2}{\rm
Im}\,k_{2}\big{(}(k-k_{*})n-2\epsilon\big{)}=\\\
&\sum\limits_{k=-\infty}^{\infty}\Big{(}X^{-3,2}_{k}(e)\Big{)}^{2}\left(\frac{h_{1}\tau_{1}\big{(}(k-k_{*})n-2\epsilon\big{)}}{1+\tau_{1}^{2}\big{(}(k-k_{*})n-2\epsilon\big{)}^{2}}+\frac{(1-h_{1})\tau_{2}\big{(}(k-k_{*})n-2\epsilon\big{)}}{1+\tau_{2}^{2}\big{(}(k-k_{*})n-2\epsilon\big{)}^{2}}\right).\end{split}$

For $k\neq k_{*}$, where $|\pm n-2\epsilon|>\text{constant}>0$, and
$\tau_{2}\gg 1$,

$\left|\frac{\tau_{2}\big{(}(k-k_{*})n-2\epsilon\big{)}}{1+\tau_{2}^{2}\big{(}(k-k_{*})n-2\epsilon\big{)}^{2}}\right|<\left|\frac{1}{\tau_{2}\big{(}(k-k_{*})n-2\epsilon\big{)}}\right|\approx
0\,,$

and the equation for the torque-free curve can be approximately written as

$\sum\limits_{k=-\infty}^{\infty}\frac{\Big{(}X^{-3,2}_{k}(e)\Big{)}^{2}}{\Big{(}X^{-3,2}_{k_{*}}(e)\Big{)}^{2}}\frac{\tau_{1}\big{(}(k-k_{*})n-2\epsilon\big{)}}{1+\tau_{1}^{2}\big{(}(k-k_{*})n-2\epsilon\big{)}^{2}}=\frac{1-h_{1}}{h_{1}}\frac{2\tau_{2}\epsilon}{1+(2\tau_{2}\epsilon)^{2}}\,,$

where, according to approximation (2.25),
$e=\sqrt{1-\left(\frac{\ell^{3}_{\scriptscriptstyle
T}}{a^{3}_{1}}n\right)^{2/3}}$. It is convenient to rewrite this equation
using the variables $\nu:=\frac{\ell^{3}_{\scriptscriptstyle
T}}{a^{3}_{1}}\,n$, ${\rm w}:=\frac{\ell^{3}_{\scriptscriptstyle
T}}{a^{3}_{1}}\,\omega$, $T_{1}:=\frac{a^{3}_{1}}{\ell^{3}_{\scriptscriptstyle
T}}\tau_{1}$, and $x=2\epsilon\tau_{2}$:

$\sum\limits_{k=-\infty}^{\infty}\frac{\Big{(}X^{-3,2}_{k}(e)\Big{)}^{2}}{\Big{(}X^{-3,2}_{k_{*}}(e)\Big{)}^{2}}\frac{\big{(}(k-k_{*})T_{1}\nu-\frac{\tau_{1}}{\tau_{2}}x\big{)}}{1+\big{(}(k-k_{*})T_{1}\nu-\frac{\tau_{1}}{\tau_{2}}x\big{)}^{2}}=\frac{1-h_{1}}{h_{1}}\frac{x}{1+x^{2}}\,,$
(3.69)

where $e=\sqrt{1-\nu^{2/3}}$.

We will assume that $0.65<\nu\leq 1$. There are two reasons for this
assumption. Firstly, we computed Hansen’s coefficients up to the 20th order in
$e$, which imposes the limitation $\nu>0.65$. The second reason is that we
will assume $\frac{\tau_{1}}{\tau_{2}}\ll T_{1}\nu$. This inequality, combined
with $|\epsilon|\tau_{2}<1$ and $|x|<2$, implies that for $k\neq k_{*}$,
$(k-k_{*})T_{1}\nu-\frac{\tau_{1}}{\tau_{2}}x\approx(k-k_{*})T_{1}\nu$.
Therefore, equation (3.69) can be approximated by

$\sum\limits_{k\neq
k_{*}}\frac{\Big{(}X^{-3,2}_{k}(e)\Big{)}^{2}}{\Big{(}X^{-3,2}_{k_{*}}(e)\Big{)}^{2}}\frac{\big{(}(k-k_{*})T_{1}\nu\big{)}}{1+\big{(}(k-k_{*})T_{1}\nu\big{)}^{2}}=\frac{1-h_{1}}{h_{1}}\frac{x}{1+x^{2}}\,,$
(3.70)

where the term
$\frac{\frac{\tau_{1}}{\tau_{2}}x}{1+\big{(}\frac{\tau_{1}}{\tau_{2}}x\big{)}^{2}}$,
corresponding to $k=k_{*}$, was neglected because
$\left|\frac{\tau_{1}}{\tau_{2}}x\right|$ is small.

The function on the right-hand side of equation (3.70) depends solely on $x$
and, except for the factor $\frac{1-h_{1}}{h_{1}}$, it is plotted in Figure
11. The function on the left-hand side of equation (3.70) depends only on
$\nu$ and $T_{1}$. For $k_{*}=4$, which corresponds to the resonance
$2n-\omega=0$, this graph with $T_{1}=0.512$ is plotted in Figure 12.

Figure 11: The function on the right-hand side of equation (3.70). Figure 12:
Graph (solid blue curve) of the function on the left-hand side of equation
(3.70) for $T_{1}=0.512$, which corresponds to $\tau_{1}n_{mer}=0.48$. The
horizontal line (brown-dashed) corresponds to
$-4.045=-\frac{1-h_{1}}{2h_{1}}$, $h_{1}=0.11$.

The function on the right-hand side of equation (3.70) has its minimum equal
to $-\frac{1-h_{1}}{2h_{1}}$ and its maximum equal to
$\frac{1-h_{1}}{2h_{1}}$. For a given value of $h_{1}\in(0,1)$, the maximum
value of $\nu$ on the torque-free curve is determined by the intersection of
the horizontal line $-\frac{1-h_{1}}{2h_{1}}$ and the graph of the function on
the left-hand side of equation (3.70), as depicted in Figure 12. This result
should be accurate in the limit as $\tau_{2}\to+\infty$ and is likely
reasonable for $\tau_{2}/\tau_{1}\gg 1$.

Using Figure 12, we estimated that on the torque-free curve, close to the
resonance $2n-\omega=0$, and for $T_{1}=0.512$ (which corresponds to
$\tau_{1}n_{\text{mer}}=0.48$), $\nu$ attains a maximum value of $0.927077$.
This value is close to $0.9267$, which is the value for the corresponding
point in Figure 6, where $\frac{\tau_{2}}{\tau_{1}}=2900$.

## Acknowledgements

L. R. S. is supported in part by FAPEMIG (Fundação de Amparo à Pesquisa no
Estado de Minas Gerais) under Grants No. RED-00133-21 and APQ-02153-23.

## References

  * Alexander [1973] ME Alexander.  The weak friction approximation and tidal evolution in close binary systems.  _Astrophysics and Space Science_ , 23:459–510, 1973.
  * Baland et al. [2017] Rose-Marie Baland, Marie Yseboodt, Attilio Rivoldini, and Tim Van Hoolst.  Obliquity of Mercury: Influence of the precession of the pericenter and of tides.  _Icarus_ , 291:136–159, 2017.
  * Bland [1960] DR Bland.  _Linear viscoelasticity_.  Pergamon Press, Oxford, 1960.
  * Boué and Efroimsky [2019] Gwenaël Boué and Michael Efroimsky.  Tidal evolution of the Keplerian elements.  _Celestial Mechanics and Dynamical Astronomy_ , 131:1–46, 2019.
  * Boué et al. [2016] Gwenaël Boué, Alexandre CM Correia, and Jacques Laskar.  Complete spin and orbital evolution of close-in bodies using a Maxwell viscoelastic rheology.  _Celestial Mechanics and Dynamical Astronomy_ , 126(1-3):31–60, 2016.
  * Consorzi et al. [2023] Anastasia Consorzi, Daniele Melini, and Giorgio Spada.  Relation between the moment of inertia and the $k_{2}$ Love number of fluid extra-solar planets.  _Astronomy & Astrophysics_, 676:A21, 2023.
  * Correia et al. [2018] ACM Correia, C Ragazzo, and LS Ruiz.  The effects of deformation inertia (kinetic energy) in the orbital and spin evolution of close-in bodies.  _Celestial Mechanics and Dynamical Astronomy_ , 130(8):51, 2018.
  * Correia and Valente [2022] Alexandre CM Correia and Ema FS Valente.  Tidal evolution for any rheological model using a vectorial approach expressed in Hansen coefficients.  _Celestial Mechanics and Dynamical Astronomy_ , 134(3):24, 2022.
  * Correia et al. [2014] Alexandre CM Correia, Gwenaël Boué, Jacques Laskar, and Adrián Rodríguez.  Deformation and tidal evolution of close-in planets and satellites using a Maxwell viscoelastic rheology.  _Astronomy & Astrophysics_, 571:A50, 2014.
  * Efroimsky [2012] Michael Efroimsky.  Bodily tides near spin–orbit resonances.  _Celestial Mechanics and Dynamical Astronomy_ , 112(3):283–330, 2012.
  * Ferraz-Mello [2013] Sylvio Ferraz-Mello.  Tidal synchronization of close-in satellites and exoplanets. a rheophysical approach.  _Celestial Mechanics and Dynamical Astronomy_ , 116(2):109–140, 2013.
  * Ferraz-Mello [2015] Sylvio Ferraz-Mello.  Tidal synchronization of close-in satellites and exoplanets: II. Spin dynamics and extension to mercury and exoplanet host stars.  _Celestial Mechanics and Dynamical Astronomy_ , 122:359–389, 2015.
  * Ferraz-Mello [2019] Sylvio Ferraz-Mello.  Planetary tides: theories.  _Satellite Dynamics and Space Missions_ , pages 1–50, 2019.
  * Ferraz-Mello [2021] Sylvio Ferraz-Mello.  On tides and exoplanets.  _Proceedings of the International Astronomical Union_ , 15(S364):20–30, 2021.
  * Ferraz-Mello et al. [2020] Sylvio Ferraz-Mello, Cristian Beaugé, Hugo A Folonier, and Gabriel O Gomes.  Tidal friction in satellites and planets. the new version of the creep tide theory.  _The European Physical Journal Special Topics_ , 229:1441–1462, 2020.
  * Fienga et al. [2019] Agnès Fienga, Pierre Deram, V Viswanathan, A Ruscio, L Bernus, D Durante, M Gastineau, and J Laskar.  _INPOP19a planetary ephemerides_.  PhD thesis, IMCCE, 2019.
  * Folonier et al. [2018] HA Folonier, S Ferraz-Mello, and E Andrade-Ines.  Tidal synchronization of close-in satellites and exoplanets. III. tidal dissipation revisited and application to Enceladus.  _Celestial Mechanics and Dynamical Astronomy_ , 130(12):78, 2018.
  * Genova et al. [2019] Antonio Genova, Sander Goossens, Erwan Mazarico, Frank G Lemoine, Gregory A Neumann, Weijia Kuang, Terence J Sabaka, Steven A Hauck, David E Smith, Sean C Solomon, et al.  Geodetic evidence that Mercury has a solid inner core.  _Geophysical Research Letters_ , 46(7):3625–3633, 2019.
  * Gevorgyan [2021] Yeva Gevorgyan.  Homogeneous model for the TRAPPIST-1e planet with an icy layer.  _Astronomy & Astrophysics_, 650:A141, June 2021.  doi: 10.1051/0004-6361/202140736.
  * Gevorgyan et al. [2020] Yeva Gevorgyan, Gwenaël Boué, Clodoaldo Ragazzo, Lucas S. Ruiz, and Alexandre C.M. Correia.  Andrade rheology in time-domain. Application to Enceladus’ dissipation of energy due to forced libration.  _Icarus_ , 343:113610, 2020.
  * Gevorgyan et al. [2023] Yeva Gevorgyan, Isamu Matsuyama, and Clodoaldo Ragazzo.  Equivalence between simple multilayered and homogeneous laboratory-based rheological models in planetary science.  _Monthly Notices of the Royal Astronomical Society_ , 523(2):1822–1831, 2023.
  * Goldreich [1966] Peter Goldreich.  Final spin states of planets and satellites.  _The Astronomical Journal_ , 71:1, 1966.
  * Gomes et al. [2019] Gabriel O Gomes, Hugo A Folonier, and Sylvio Ferraz-Mello.  Rotation and figure evolution in the creep tide theory: a new approach and application to Mercury.  _Celestial Mechanics and Dynamical Astronomy_ , 131(12):56, 2019.
  * Goossens et al. [2022] Sander Goossens, Joe P Renaud, Wade G Henning, Erwan Mazarico, Stefano Bertone, and Antonio Genova.  Evaluation of recent measurements of Mercury’s moments of inertia and tides using a comprehensive Markov Chain Monte Carlo method.  _The Planetary Science Journal_ , 3(2):37, 2022\.
  * Hut [1981] P Hut.  Tidal evolution in close binary systems.  _Astronomy and Astrophysics, vol. 99, no. 1, June 1981, p. 126-140._ , 99:126–140, 1981.
  * Kaula [1964] William M Kaula.  Tidal dissipation by solid friction and the resulting orbital evolution.  _Reviews of geophysics_ , 2(4):661–685, 1964\.
  * Makarov [2012] Valeri V Makarov.  Conditions of passage and entrapment of terrestrial planets in spin-orbit resonances.  _The Astrophysical Journal_ , 752(1):73, 2012\.
  * Makarov and Efroimsky [2013] Valeri V Makarov and Michael Efroimsky.  No pseudosynchronous rotation for terrestrial planets and moons.  _The Astrophysical Journal_ , 764(1):27, 2013\.
  * Margot et al. [2018] Jean-Luc Margot, SA Hauck II, E Mazarico, Stanton J Peale, and Sebastiano Padovan.  Mercury’s internal structure.  _Mercury-The view after MESSENGER_ , pages 85–113, 2018.
  * Matsuyama and Nimmo [2009] I Matsuyama and F Nimmo.  Gravity and tectonic patterns of Mercury: Effect of tidal deformation, spin-orbit resonance, nonzero eccentricity, despinning, and reorientation.  _Journal of Geophysical Research: Planets_ , 114(E1), 2009\.
  * Matsuyama et al. [2016] Isamu Matsuyama, Francis Nimmo, James T Keane, Ngai H Chan, G Jeffrey Taylor, Mark A Wieczorek, Walter S Kiefer, and James G Williams.  GRAIL, LLR, and LOLA constraints on the interior structure of the Moon.  _Geophysical Research Letters_ , 43(16):8365–8375, 2016.
  * Mignard [1979] F Mignard.  The evolution of the lunar orbit revisited. I.  _The Moon and the planets_ , 20(3):301–315, 1979\.
  * Noyelles et al. [2014] Benoit Noyelles, Julien Frouard, Valeri V Makarov, and Michael Efroimsky.  Spin–orbit evolution of Mercury revisited.  _Icarus_ , 241:26–44, 2014.
  * Ogilvie [2014] Gordon I Ogilvie.  Tidal dissipation in stars and giant planets.  _Annual Review of Astronomy and Astrophysics_ , 52:171–210, 2014.
  * Ragazzo and Ruiz [2024] C Ragazzo and L S Ruiz.  Spin-orbit synchronization and singular perturbation theory.  _To appear_ , 2024.
  * Ragazzo and Ruiz [2015] C Ragazzo and LS Ruiz.  Dynamics of an isolated, viscoelastic, self-gravitating body.  _Celestial Mechanics and Dynamical Astronomy_ , 122(4):303–332, 2015.
  * Ragazzo [2020] Clodoaldo Ragazzo.  The theory of figures of Clairaut with focus on the gravitational modulus: inequalities and an improvement in the Darwin–Radau equation.  _São Paulo Journal of Mathematical Sciences_ , 14:1–48, 2020.
  * Ragazzo and Ruiz [2017] Clodoaldo Ragazzo and LS Ruiz.  Viscoelastic tides: models for use in Celestial Mechanics.  _Celestial Mechanics and Dynamical Astronomy_ , 128(1):19–59, 2017.
  * Ragazzo et al. [2022] Clodoaldo Ragazzo, Gwenaël Boué, Yeva Gevorgyan, and Lucas S Ruiz.  Librations of a body composed of a deformable mantle and a fluid core.  _Celestial Mechanics and Dynamical Astronomy_ , 134(2):10, 2022.
  * Rochester and Smylie [1974] MG Rochester and DE Smylie.  On changes in the trace of the Earth’s inertia tensor.  _Journal of Geophysical Research_ , 79(32):4948–4951, 1974.
  * Sabadini et al. [2016] Roberto Sabadini, Bert Vermeersen, and Gabriele Cambiotti.  _Global dynamics of the Earth_.  Springer, 2016.
  * Singer [1968] SF Singer.  The origin of the Moon and geophysical consequences.  _Geophysical Journal International_ , 15(1-2):205–226, 1968.
  * Steinbrügge et al. [2018] G Steinbrügge, S Padovan, H Hussmann, T Steinke, A Stark, and J Oberst.  Viscoelastic Tides of Mercury and the Determination of its Inner Core Size.  _Journal of Geophysical Research: Planets_ , 123(10):2760–2772, 2018.