11institutetext: Observatorio Astronómico de Córdoba, Universidad Nacional de
Córdoba, Laprida 854, Córdoba X5000BGR, Argentina
11email<EMAIL_ADDRESS>22institutetext: CONICET, Instituto de
Astronomía Teórica y Experimental, Laprida 854, Córdoba X5000BGR, Argentina
33institutetext: Instituto de Astronomia Geofísica e Ciências Atmosféricas,
Universidade de São Paulo, SP 05508-090, Brazil

# Creep Tide Model for the 3-Body Problem

The rotational evolution of a circumbinary planet

F. A. Zoppetti 1122 H. Folonier 3344 A. M. Leiva 11 C. Beaugé 1122

We present a tidal model for treating the rotational evolution in the general
three-body problem with arbitrary viscosities, in which all the masses are
considered to be extended and all the tidal interactions between pairs are
taken into account. Based on the creep tide theory, we present the set of
differential equations that describes the rotational evolution of each body,
in a formalism that is easily extensible to the N tidally-interacting body
problem. We apply our model to the case of a circumbinary planet and use a
Kepler-38 like binary system as a working example. We find that, in this low
planetary eccentricity case, the most likely final stationary rotation state
is the 1:1 spin-orbit resonance, considering an arbitrary planetary viscosity
inside the estimated range for the solar system planets. The timescales for
reaching the equilibrium state are expected to be $\sim$ Myrs for stiff bodies
but can be longer than the age of the system for planets with a large gaseous
component. We derive analytical expressions for the mean rotational stationary
state, based on high-order power series of the semimajor axes ratio
$a_{1}/a_{2}$ and low-order expansions of the eccentricities. These are found
to reproduce very accurately the mean behaviour of the low-eccentric numerical
integrations for arbitrary planetary relaxation factors, and up to
$a_{1}/a_{2}\sim 0.4$. Our analytical model is used to predict the stationary
rotation of the Kepler circumbinary planets and find that most of them are
probably rotating in a sub-synchronous state, although the synchrony shift is
much less important than the one estimated in Zoppetti et al. (2019, 2020). We
present a comparison of our results with those obtained with the Constant Time
Lag and find that, unlike what we assumed in our previous works, the cross
torques have a non-negligible net secular contribution, and must be taken into
account when computing the tides over each body in an N-extended-body system
from an arbitrary reference frame. These torques are naturally taken into
account in the creep theory. In addition to this, the latter formalism
considers more realistic rheology that proved to reduce to the Constant Time
Lag model in the gaseous limit and also allows to study several additional
relevant physical phenomena.

###### Key Words.:

planet-star interactions – planets and satellites: dynamical evolution and
stability – methods: analytical

##  1 Introduction

It is well known that the rotational evolution of extended bodies in compact
systems is strongly affected by tides. In the two-body problem (2BP), tides
tend to drive the spins to a stationary state characterized by a spin-orbit
resonance. In particular, tidal models following the Darwin theory predict the
existence of a stationary solution synchronous with the mean motion, when the
two bodies move in circular orbits, or super-synchronous otherwise (e.g.
Ferraz-Mello et al., 2008). However, this prediction has not been confirmed by
the observations of satellites rotation rates, as in the case of Titan (e.g.
Meriggiola 2012), Europa (e.g. Hurford et al., 2007) and the moon itself.

Several attempts were made to explain these differences, either considering an
extra torque caused by the existence of a permanent deformation, or by
modifying the tidal theories. In particular, Ferraz-Mello (2013, 2015)
proposed a theory of tidal interactions based on an approximate solution of
the Navier–Stokes equation, where the capture in spin–orbit resonances can be
explained without the assumption of an extra torque. This result was also
reported by Makarov & Efroimsky (2013) and confirmed by Correia et al. (2014)
with a model based on a Maxwell viscoelastic rheology.

Regarding the general three-body problem (3BP), a model that takes into
account all tidal interactions between pairs and valid for bodies with
arbitrary viscosities is still missing. Instead, this problem is usually
addressed by neglecting some of the tidal interactions. For example, in the
circumbinary problem in which the planetary mass is expected to be much
smaller than the stellar masses, the planetary tides induced on the stellar
components are commonly omitted (e.g. Correia et al., 2016). Although the
appropriateness of these approximations is not under debate, building a
general model valid for any mass ratio between pairs, allows extending the
results to problems with different hierarchies, such as triple star systems
and exomoons host planets close to the central star.

In Zoppetti et al. (2019) we presented a self-consistent tidal model for the
planar and non-oblique 3BP, using a weak friction Constant Time Lag (CTL)
model, based on the expressions of the tidal forces and torques of the
classical work of Mignard (1979). We applied the model to the circumbinary
(CB) case and found that for systems like those discovered by the Kepler
mission, the planet acquired stationary spin rates on timescales of $\sim$
Myrs (depending on the dissipation factor $Q_{2}$) and, curiously, in sub-
synchronous solutions (Zoppetti et al., 2020). The sub-synchronous of the
planetary spin was found to be the result of the tidal interaction with the
secondary central star and proportional to its mass, while the planetary
eccentricity tended to generate super-synchronous states, as in the 2BP.
However, the goal of building a general tidal model in the 3BP has not yet
been fulfilled since the weak-friction regime of the CTL models has shown to
be only valid for bodies with a large gaseous component (Efroimsky, 2012,
2015).

In this work, we present and discuss a self-consistent model for treating the
tides in the 3BP, in the framework of the Newtonian creep tidal model Ferraz-
Mello (2013). We focus our attention on the rotational evolution of the bodies
and apply the results to the case of a CB planet. In particular, to allow for
comparison with previous results, we employ a Kepler-38 like system (Orosz et
al., 2012). However, a wider range of system parameters as well as different
initial orbital elements are also explored.

This paper is organized as follows. In Section 2 we present the tidal model by
computing the equilibrium figures of the bodies and deriving differential
equations from the creep equation to model the real shapes and orientations
for each one. Once the real figures are obtained, the rotational evolution is
calculated from the reaction torques on the two companions. Section 3 presents
a series of numerical integrations of the shape, orientation, and spin
evolution equations for a CB planet in a Kepler-38 system, considering
different planetary internal viscosities. In Section 4, we construct
analytical expressions for the mean values of the rotational stationary state
of the CB planet from a Jacobian reference frame and compare these predictions
with numerical solutions. In addition to this, we estimate the equilibrium
rotation of the Kepler CB planets. In Section 5, we compare our results with
those obtained in Zoppetti et al. (2019, 2020), and present an improved
formalism that takes into account the cross torques in the CTL model. Finally,
Section 6 summarises our main results and discusses their implications.

##  2 The Creep Tidal Model for the 3BP

Let us consider the extended three-body problem with masses $m_{i}$, $m_{j}$
and $m_{k}$ and respective mean radius $\mathcal{R}_{i}$, $\mathcal{R}_{j}$
and $\mathcal{R}_{k}$. In addition, let us assume that the orbital motion of
the bodies lies in the same plane and the spin vectors are perpendicular to
it. From an arbitrary inertial frame, let ${\bf R}_{i}$, ${\bf R}_{j}$ and
${\bf R}_{k}$ be the respective position vectors of $m_{i}$, $m_{j}$ and
$m_{k}$. We denote the relative position of $m_{j}$ to $m_{i}$ by
${\bf\Delta}_{ji}={\bf R}_{j}-{\bf R}_{i}$, while the angle subtended between
the position of $m_{j}$ and the reference axis, as seen from $m_{i}$, is
called $\varphi_{ji}$ (the same for $j=k$, see Figure 1).

Figure 1: The three-body problem from an inertial frame and from a
$m_{i}$-relative frame.

Our goal is to study the spin dynamics of these three bodies under the effects
of their mutual tidal perturbations. The tides will be modeled employing the
creep theory (Ferraz-Mello, 2013, 2015) in its closed form (Folonier et al.,
2018) and extended to the case of three interacting masses.

###  2.1 Equilibrium figure

Figure 2: Left: Equilibrium figure in the 3BP. Dashed-empty ellipsoids
correspond to the static tidal deformations generated by each of the
perturbing masses, while the resulting equilibrium figure of $m_{i}$ (equation
(3)) is shown by a filled gray ellipsoid. The unit vector ${\hat{\bf
r}}^{(eq)}_{i}$ indicates the orientation of the equivalent equilibrium bulge.
Right: Real figure in the 3BP. The dark grey filled ellipsoid represents the
modelled real figure induced on $m_{i}$ when we take into account the internal
viscosity. Note that it is displaced by an angle $\delta_{i}$ respect to the
equilibrium figure and also has different flattenings.

As a consequence of the gravitational interaction on extended bodies, each
mass $m_{i}$ undergoes a deformation which is given by the superposition of
the tidal interaction with the two companions $m_{j}$ and $m_{k}$. In this
work, we also consider the rotational flattening on each body due to its own
spin, and assume that the resulting deformation due to these effects is small
enough that a model developed up to the first order of the flattenings (like
the one we will develop in this work) represents an accurate description of
the problem.

If we consider that each mass $m_{i}$ is a perfect fluid that instantly
responds to any external stress, it will immediately adopt an equilibrium
figure determined by the tides raised by the two companions (i.e. static
tides) and its own rotation. The surface equation of this equilibrium figure
seen from the center of mass, is determined by giving the distance $\rho_{i}$
of an arbitrary point on its surface with co-latitude angle $\theta_{i}$ and
longitude angle $\varphi_{i}$. It can be expressed as the sum of the
interactions between pairs as (Tisserand, 1891; Folonier et al., 2015)

$\begin{split}\rho_{i}({\theta_{i}},{\varphi_{i}})&=\mathcal{R}_{i}\bigg{[}1+\bigg{(}\epsilon_{i}^{M}+\frac{\epsilon_{ij}^{\rho}}{2}+\frac{\epsilon_{ik}^{\rho}}{2}\bigg{)}\biggl{(}\frac{1}{3}-\cos^{2}{\theta_{i}}\biggr{)}\\\
&+\frac{1}{2}\sin^{2}{\theta_{i}}\biggl{(}\epsilon^{\rho}_{ij}\cos{(2{\varphi_{i}}-2\varphi_{ji})}+\epsilon^{\rho}_{ik}\cos{(2{\varphi_{i}}-2\varphi_{ki})}\biggr{)}\bigg{]}\end{split}$
(1)

where we have chosen the same definitions for the flattening as in Folonier et
al. (2018), given by

$\displaystyle\epsilon_{ij}^{\rho}=\frac{15m_{j}\mathcal{R}_{i}^{3}}{4m_{i}\Delta^{3}_{ji}}\
\ $ , $\displaystyle\ \
\epsilon_{i}^{M}=\frac{5\mathcal{R}_{i}^{3}\Omega_{i}^{2}}{4m_{i}\mathcal{G}},$
(2)

with $\Omega_{i}$ the spin module of $m_{i}$ and $\mathcal{G}$ the
gravitational constant.

We propose that the equilibrium figure expressed in equation (1) corresponds
to an equivalent triaxial flat ellipsoid with equatorial flattening
$\epsilon_{i}^{\rho}$, polar flattening $\epsilon_{i}^{z}$ and equilibrium
orientation defined by the unit vector ${\hat{\bf
r}}^{(eq)}_{i}=(\cos{\varphi^{(eq)}_{i}},\sin{\varphi^{(eq)}_{i}})$ (see left
panel of Figure 2). The surface equation of this equivalent equilibrium figure
is simply given by

$\begin{split}\rho_{i}({\theta_{i}},{\varphi_{i}})&=\mathcal{R}_{i}\bigg{[}1+\epsilon^{z}_{i}\bigg{(}\frac{1}{3}-\cos^{2}{\theta_{i}}\bigg{)}+\frac{\epsilon^{\rho}_{i}}{2}\sin^{2}{\theta_{i}}\cos{(2\varphi_{i}-2\varphi_{i}^{eq})}\biggr{]}\end{split}$
(3)

The equivalent flattenings and orientation can be obtained in terms of the
two-body interacting parameters by equating equation (3) and equation (1).
This results in

$\displaystyle\epsilon_{i}^{z}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\bigg{(}\epsilon_{ij}^{\rho}+\epsilon_{ik}^{\rho}\bigg{)}+\epsilon_{i}^{M}$
$\displaystyle{(\epsilon_{i}^{\rho})}^{2}$ $\displaystyle=$
$\displaystyle{(\epsilon_{ij}^{\rho})}^{2}+{(\epsilon_{ik}^{\rho})}^{2}+2\epsilon_{ij}^{\rho}\epsilon_{ik}^{\rho}\cos{(2\varphi_{ji}-2\varphi_{ki})}$
(4) $\displaystyle\tan{(2\varphi_{i}^{eq})}$ $\displaystyle=$
$\displaystyle\frac{\epsilon_{ij}^{\rho}\sin{(2\varphi_{ji})}+\epsilon_{ik}^{\rho}\sin{(2\varphi_{ki})}}{\epsilon_{ij}^{\rho}\cos{(2\varphi_{ji})}+\epsilon_{ik}^{\rho}\cos{(2\varphi_{ki})}}.$

###  2.2 Creep equations

Now, we assume the more realistic case in which each body has a certain
viscosity that does not allow it to attain the equilibrium figure instantly.
The real shape of $m_{i}$, measured by the surface distance
$\zeta_{i}({\theta_{i}},{\varphi_{i}})$, is also modelled as a triaxial flat
ellipsoid but characterized by the equatorial flattening
${\mathcal{E}}_{i}^{\rho}$ and polar flattening ${\mathcal{E}}_{i}^{z}$. On
the other hand, the orientation of the real figure is displaced respect to the
equilibrium orientation by a lag $\delta_{i}$. Therefore, the real orientation
unit vector is ${\hat{\bf
r}}^{(b)}_{i}=(\cos{(\varphi^{(eq)}_{i}+\delta_{i})},\sin{(\varphi^{(eq)}_{i}+\delta_{i})})$
(see right panel in Figure 2). Thus, the surface equation of the real flat
ellipsoid is then given by

$\zeta_{i}(\theta_{i},\varphi_{i})=\mathcal{R}_{i}\bigg{[}1+{\mathcal{E}}_{i}^{z}\bigg{(}\frac{1}{3}-\cos^{2}{\theta_{i}}\bigg{)}+\frac{{\mathcal{E}}_{i}^{\rho}}{2}\sin^{2}{\theta_{i}}\cos{(2\varphi_{i}-2(\varphi_{i}^{eq}+\delta_{i}))}\bigg{]}.$
(5)

Contrary to Darwin-based tidal models, in the creep tide theory we assume that
$({\mathcal{E}}_{i}^{\rho},{\mathcal{E}}_{i}^{z},\delta_{i})$ are a-priori
unknown quantities, whose values at each instant of time are obtained by
solving the creep equation:

$\displaystyle\dot{\zeta}_{i}+\gamma_{i}\zeta_{i}=\gamma_{i}\rho_{i},$ (6)

where $\gamma_{i}$ is the so-called relaxation factor of $m_{i}$ (Ferraz-
Mello, 2013) which is inversely proportional to the body viscosity and will be
assumed as constant in this work.

A set of three differential equations for the shape and orientation of the
real figure of $m_{i}$ can be obtained by differentiating equation (5) with
respect to time, introducing its expression plus equations (3) and (5) into
equation (6) and identifying terms of equal dependence with
$(\theta_{i},\varphi_{i})$ (see Folonier et al. (2018)). Particularly, it is
convenient to express the equatorial flattening and orientation in terms of
regular variables as in Gomes et al. (2019) as

$\displaystyle\dot{X_{i}}$
$\displaystyle=-\gamma_{i}X_{i}-2Y_{i}(\Omega_{i}-\dot{\varphi_{i}}^{eq})+\gamma_{i}\epsilon_{i}^{\rho}$
$\displaystyle\dot{Y_{i}}$
$\displaystyle=-\gamma_{i}Y_{i}+2X_{i}(\Omega_{i}-\dot{\varphi_{i}}^{eq})$ (7)
$\displaystyle\dot{{\mathcal{E}}_{i}^{z}}$
$\displaystyle=\gamma_{i}({\epsilon_{i}}^{z}-{{\mathcal{E}}_{i}}^{z})$

where the regular variables are defined as

$\displaystyle X_{i}={\mathcal{E}}_{i}^{\rho}\cos(2\delta_{i})\ \ \ \ $ ;
$\displaystyle\ \ \ \ Y_{i}={\mathcal{E}}_{i}^{\rho}\sin(2\delta_{i}).$ (8)

By comparing equations (2.2) with the shape and orientation evolution
equations in the 2BP (e.g. Folonier et al., 2018; Gomes et al., 2019), we note
that in the 3BP the creep equations maintain the same form, which make the
formalism applicable to any N-body problem. Thus, the complexity in the multi-
body problem is only reduced to calculate the equilibrium figures.

###  2.3 Tidal forces and torques

At any instant of time, given the rotation frequency $\Omega_{i}$ of the body
$m_{i}$, it is possible to calculate its real shape and orientation by solving
the differential equation system $(\ref{eq:forma})$. With these quantities, it
is then possible to calculate the tidal forces and torques between each pair
of bodies.

The tidal potential due to the deformation of $m_{i}$, induced on $m_{j}$ can
be calculated as (e.g. Murray & Dermott, 1999)

$\begin{split}U_{ji}({\bf{\Delta}}_{ji})&=-\frac{{\cal
G}(B_{i}-A_{i})}{2\Delta_{ji}^{3}}\Big{(}3({{\bf{\hat{\Delta}}}_{ji}\cdot{\bf{\hat{r}}}_{i}^{b}})^{2}-1\Big{)}-\frac{{\cal
G}(C_{i}-B_{i})}{2\Delta_{ji}^{3}},\end{split}$ (9)

where $A_{i},B_{i},C_{i}$ are the principal moments of inertia of $m_{i}$,
which can be related to the ellipsoidal flattenings (see Appendix A of
Folonier et al. (2015)) as:

$(B_{i}-A_{i})\simeq\overline{\mathcal{C}}_{i}{\mathcal{E}}^{\rho}_{i}\hskip
14.22636pt;\hskip
14.22636pt(C_{i}-B_{i})\simeq\overline{\mathcal{C}}_{i}\biggl{(}{\mathcal{E}}^{z}_{i}-\frac{1}{2}{\mathcal{E}}^{\rho}_{i}\biggr{)},$
(10)

where $\overline{\mathcal{C}}_{i}=(2/5)m_{i}\mathcal{R}_{i}^{2}$ is the
principal moment of inertia of a homogeneous spherical body of mass $m_{i}$
and radius $\mathcal{R}_{i}$.

The tidal force applied on $m_{j}$ due to the deformation of $m_{i}$ can be
simply calculated as:

$\displaystyle{\bf{F}}_{ji}$ $\displaystyle=$ $\displaystyle-
m_{j}\nabla_{{\bf\Delta}_{ji}}U_{ji}$ (11) $\displaystyle=$
$\displaystyle\frac{3{\cal
G}m_{j}\overline{\mathcal{C}}_{i}}{2\Delta_{ji}^{5}}{\mathcal{E}}^{\rho}_{i}\sin{(2(\varphi^{eq}_{i}+\delta_{i})-2\varphi_{ji})}\
(\hat{\bf{k}}\times{\bf{\Delta}}_{ji})$ $\displaystyle-$
$\displaystyle\frac{3{\cal
G}m_{j}\overline{\mathcal{C}}_{i}}{2\Delta_{ji}^{5}}\Big{(}\frac{3}{2}{\mathcal{E}}^{\rho}_{i}\cos{(2(\varphi^{eq}_{i}+\delta_{i})-2\varphi_{ji})}+{\mathcal{E}}^{z}_{i}\Big{)}{\bf{\Delta}}_{ji},$

where we have decomposed $\hat{\bf{r}}_{i}^{b}$ in terms of the radial unit
vector $\hat{\bf{\Delta}}_{ji}$ and the azimutal unit vector
$(\hat{\bf{k}}\times\hat{\bf{\Delta}}_{ji})$ ($\hat{\bf{k}}$ is the unit
vector paralell to the spin vectors).

The torque induced on $m_{j}$ due to the deformation of $m_{i}$ can be simply
calculated as

$\begin{split}{\bf{T}}_{ji}&={\bf{\Delta}}_{ji}\times{\bf{F}}_{ji}\\\
&=\frac{3{\cal
G}m_{j}\overline{\mathcal{C}}_{i}}{2\Delta_{ji}^{3}}{\mathcal{E}}^{\rho}_{i}\sin{(2(\varphi^{eq}_{i}+\delta_{i})-2\varphi_{ji})}\
\hat{\bf{k}},\end{split}$ (12)

while a completely analog expression can be found for the torque applied on
$m_{k}$ due to the deformation of $m_{i}$ by exchanging the indexes.

Now, we assume that the reaction torque due to the deformation on $m_{i}$ is
the only responsable in changing its own rotational velocity. This is a
natural assumption that satisfies the condition that the variation of the
rotation speed of $m_{i}$ should be cancelled in the limit $\mathcal{R}_{i}\to
0$ (Zoppetti et al., 2019). Thus, the rotational evolution of $m_{i}$ is given
by

$\displaystyle C_{i}\dot{{\bf{\Omega}}}_{i}+\dot{C}_{i}{{\bf{\Omega}}}_{i}$
$\displaystyle=$ $\displaystyle-\Big{(}{{\bf{T}}}_{ji}+{{\bf{T}}}_{ki}\Big{)}$
(13) $\displaystyle=$ $\displaystyle-\frac{3{\cal
G}\overline{\mathcal{C}}_{i}{\mathcal{E}}^{\rho}_{i}}{2}\bigg{[}\frac{m_{j}}{\Delta_{ji}^{3}}\sin{(2(\varphi^{eq}_{i}+\delta_{i})-2\varphi_{ji})}$
$\displaystyle+$
$\displaystyle\frac{m_{k}}{\Delta_{ki}^{3}}\sin{(2(\varphi^{eq}_{i}+\delta_{i})-2\varphi_{ki})}\
\bigg{]}\hat{\bf{k}}$

Neglecting the term corresponding to the variation of the polar moment of
inertia (i.e. $\dot{C}_{i}\sim 0$), using equations (2.1) and after some
algebra, it is possibly to find that the rotational evolution of the body
$m_{i}$ is governed by

$\displaystyle\dot{\Omega}_{i}$ $\displaystyle=-\frac{2{\cal
G}m_{i}}{5\mathcal{R}_{i}^{3}}\epsilon_{i}^{(\rho)}Y_{i}$ (14)

The set of differential equations formed by equations (2.2) and (14) described
the dynamical evolution of the system on short timescales and will be target
of this work.

##  3 Numerical simulations

Circumbinary planetary systems are particularly interesting scenarios where we
can study the effects of tides on the dynamical evolution, in the framework of
the 3BP. In this context, the stellar tidal evolution is expected to be very
little affected by the presence of the CB planet and mainly dominated by the
2-body interaction between the binary components (Correia et al., 2016;
Zoppetti et al., 2019). Since the rotational evolution due to tides in a
binary system has been extensively studied (e.g. Mignard, 1979; Hut, 1981;
Efroimsky, 2012), in this paper we will focus on the evolution of the planet.

Using a CTL model based on the tidal forces and torques of Mignard (1979), in
Zoppetti et al. (2019, 2020) we found that the characteristic timescales of
the rotational evolution of any planet in a Kepler-like circumbinary system
are much faster than the characteristic timescales of the orbital evolution of
the system. As a result, in this framework, it is expected that while the
shape, orientation and rotation of each body reach stationary solution, only
periodic oscillations be observed in the orbital elements of the bodies,
except for the fast angles. For this reason, in this work we take advantage of
this adiabatic nature of the problem and solve the set of differential
equation associated with the evolution of the shape of each body, its
orientation and its rotation (equations (2.2) and (14)), assuming fixed values
for all the orbital elements except the true anomalies. The orbital evolution
of the bodies on much longer timescales will be addressed in a forthcoming
work.

With this motivation, we begin first by analyzing the results of the numerical
integrations of the equations (2.2) and (14) in a 3BP consisting of a single
planet around a binary star. To compare with the results of our previous works
(Zoppetti et al., 2018, 2019), we choose a Kepler-38 like system as a test
case. Nominal values for system parameters and initial orbital elements are
detailed in Table 1. Stellar masses and radii were taken close to the observed
by Orosz et al. (2012), while the value of $m_{2}$ was estimated from a semi-
empirical mass-radius fit (Mills & Mazeh, 2017).

Table 1: Initial conditions of our reference numerical simulation, representing a Kepler-38 like system (Orosz et al., 2012). Orbital elements are given in a Jacobi reference frame. The parameters highlighted with an asterisk were varied in different simulations as is indicated in the text. body | $m_{0}$ | $m_{1}$ | $m_{2}$
---|---|---|---
mass | $0.95\,m_{\odot}$ |  $0.25\,m_{\odot}$(*) | $10\,m_{\oplus}$
radius | $0.84\,\mathcal{R}_{\odot}$ | $0.27\,\mathcal{R}_{\odot}$ | $4.35\,\mathcal{R}_{\oplus}$
$a_{i}$ [AU] |  | $0.15$ |  $0.45$(*)
$e_{i}$ |  | $0.15$ |  $0.05$(*)
$\Omega_{i}$ | $5\,n_{1}$ | $5\,n_{1}$ | $5\,n_{2}$
Figure 3: Rotational evolution due to tides in a Kepler-38 like CB planet (see
Table 1), considering different values of planetary relaxation factor
$\gamma_{2}$ with different colors (see legend in the upper left panel). The
evolution was calculated by direct numerical integration of equations (2.2)
and (14). The evolution of the planetary spin rate $\Omega_{2}$ is shown in
upper left panel (zoomed around $\Omega_{2}\sim n_{2}$ in a secondary plot),
the lag-angle $\delta_{2}$ in upper right panel, the equatorial flattening
$\mathcal{E}^{\rho}_{2}$ in the lower left panel and the polar flattening
$\mathcal{E}^{z}_{2}$ in the lower right panel. The dashed black line in the
upper left panel shows the stationary solution predicted by Zoppetti et al.
(2019) with a CTL-based model.

The first obstacle we encountered when doing the numerical integration was the
virtual impossibility of solving the equations system numerically for bodies
with large relaxation factors ($\gamma_{2}>>n_{1},n_{2}$), typically
associated with bodies with a large gaseous component (Ferraz-Mello, 2013).
This limitation comes from the fact that, for large $\gamma$-values, the shape
and orientation of the deformed body quickly adapt to the change in position
of the deforming bodies, so the integrator must greatly reduce the time-step
to pursue it and, finally, integration over millions of years timescales
becomes computationally impracticable (e.g. Lambert, 1991). However, as we
will see, the rotational evolution timescales can be extrapolated from the
simulations of more viscous bodies in the same regime, and the synchronous
attractor from the analytical solution calculated in Section 4.

In Figure 3 we show the rotational evolution of a CB planet from a Jacobi
reference frame, obtained by direct numerical integration of equations (2.2)
and (14). All initial conditions and system parameters were taken equal to the
nominal values summarized in Table 1, while the angles were set equal to zero.
As we mentioned above, the orbital elements are considered fixed, with the
exception of the mean anomalies which vary linearly with the time.

Many interesting features can be noticed in Figure 3. We begin first by
identifying, as in the 2BP, two different regimes whose boundary value is
close to the main dynamical frequency of the problem $n_{2}$: the ”stiff
regime” for $\gamma_{2}<n_{2}$ and the ”gaseous regime” for
$\gamma_{2}>n_{2}$. In the stiff regime, the rotational evolution acquires its
stationary value faster the higher the relaxation factor $\gamma_{2}$, while
the opposite occurs in the gaseous limit: if we decrease the planetary
viscosity (increase $\gamma_{2}$) the attractor is reached on longer
timescales. This bimodal character of the energy dissipation has already been
observed previously in the 2BP (Ferraz-Mello, 2013; Correia et al., 2014): for
relaxation parameters far from the maximum, the energy dissipation decreases
proportional to $\gamma_{2}$ in the gaseous regime and decreases inversely
proportional to $\gamma_{2}$ in the stiff bodies limit.

In the stiff regime, it is possible to notice from the upper left panel of
Figure 3 that for bodies with very small relaxation factors
$\gamma_{2}/n_{2}\sim 0.01$, ( i.e. $\gamma_{2}\simeq 7\times 10^{-9}$ Hz),
the planetary spin captures may occur in an attractor different than 1:1 spin-
orbit resonance, in our testing case 3:2 with an associated circularization of
the lag-angle $\delta_{2}$. This result is expected from the 2BP experience,
where for very low $\gamma_{2}$-values, the probability of capture in this
type of resonance increases considerably (Correia et al., 2014; Gomes et al.,
2019). However, according to our numerical simulations, captures in resonances
different than the 1:1 seems to occur for bodies with satellite-like
relaxation factors like the moon, which are expected to be approximately two
orders of magnitude lower than those expected for rocky earth-type planets
(Ferraz-Mello, 2013), or for bodies with high initial eccentricities, a case
that we will not consider in this work. This particular stationary state also
has a counterpart in the flattenings $\mathcal{E}^{\rho}_{2}$ and
$\mathcal{E}^{z}_{2}$, which oscillate around values significantly different
than in the case of the capture in the 1:1 spin-orbit resonance (see lower
panels in Figure 3).

The oscillation of the lag-angle around $\delta_{2}\sim 45^{\circ}$ of free
rotating bodies in the stiff regime can be observed in the upper right panel
of Figure 3 (blue and orange dots). This result also has been previously
reported in the framework of the 2BP (Gomes et al., 2019) and is mainly due to
the lack of the elastic tide component inside the creep model (Ferraz-Mello,
2015). Finally, we notice that in this limit, the oscillation of $\delta_{2}$
(and also of $\Omega_{2}$, although it can not be distinguished in the upper
left panels) around the stationary values are of much greater amplitudes
respect to those in the gaseous regime, while the completely opposite
situation occurs for $\mathcal{E}^{\rho}_{2}$ and $\mathcal{E}^{z}_{2}$. As
noticed by Folonier et al. (2018) in the 2BP, this periodic oscillations of
large amplitudes must be taken into account when constructing an analytical
approach (see Section 4).

Regarding the gaseous regime, despite considering bodies of different
viscosities (i.e. relaxation factors $\gamma_{2}$), we observe similar
behaviors: stationary spins are located very close to $n_{2}$ and the lag-
angles $\delta_{2}$ oscillates close to zero, with low oscillation amplitudes.
If we omit the case of $\gamma_{2}/n_{2}=1$ which is a behavior change value,
we can observe in the secondary plot of the upper left panel that the
equilibrium spins become smaller as we increase $\gamma_{2}$ and turn out to
be sub-synchronous for these parameter system in the gaseous regime (see red
and magenta curves). However, despite the fact that the rotational equilibrium
state resulting from the numerical integrations appears to be sub-synchronous
in the gaseous limit, this solution is far from the one we found in Zoppetti
et al. (2019) (dashed black line in the upper left panel). The details of the
1:1 spin-orbit mean solution will be described in Section 4, while the reasons
for the discrepancy with Zoppetti et al. (2019) will be discussed in Section
5.

##  4 An analytical solution for the mean rotational stationary state

As we have discussed in the previous Section 3, according to our numerical
simulations, the rotational evolution due to tides of a CB planet in a low-
eccentric orbit with $\gamma_{2}\gtrsim 0.1n_{2}$ (approximately the expected
range for a planet-like body, according to Ferraz-Mello (2013)) tends to be
captured in the 1:1 spin-orbit resonance.

The timescales for the capture depend on the planetary viscosities. While for
bodies in the stiff regime are expected to be $\lesssim$ Myrs, in the gaseous
regime can be of the order of some Gyrs, considering $\gamma_{2}$-value
similar to the one estimated for Neptune (Ferraz-Mello, 2013) and
extrapolating the behavior observed in the Figure 3 by assuming a linear
dependence of the decay times with the relaxation factors. Although for giant
planets these timescales may be even larger than the age of the systems, the
enormous uncertainty that exists about the internal structure of the
exoplanets, and in particular the CB Kepler planets, prohibits a reliable
estimate of the planetary relaxation factor. In addition to this, the central
binary could have undergone an important change in both the physical
parameters and the orbital elements, which may have accelerated the tidal
evolution processes (Zoppetti et al., 2018). In summary, according to our
model, a CB planet in a Kepler-38 like system with a large rocky component is
expected to have reached its stationary spin and this is also possible for the
case of gaseous CB planets with small relaxation factors $\gamma_{2}$ .

From a practical point of view, having an analytical approach for the
rotational stationary state will allow us, when studying long-term orbital
evolution of the system, to avoid direct numerical integration of the
equations (2.2) and (14), which are particularly difficult to solve for bodies
with large relaxation factors (see Section 3) and, instead, to work with a
semi-analytical approach in which we introduce a recipe for the equilibrium
shape, orientation and spin.

With these motivations, in this Section we develop an analytical approximation
for the equations system (2.2) and (14) around the pseudo-synchronous
stationary state, in this particular scenario where there are two central
bodies (both stellar components of the binary) exerting tides onto a planet.
Mathematically speaking, our goal is to find the stationary solution for the
system formed by $(X_{2},Y_{2},{\mathcal{E}}^{z}_{2})$ and $\Omega_{2}$,
employing a power series truncated in $\alpha=a_{1}/a_{2}$, $e_{1}$ and
$e_{2}$.

To do this, we first adopt a Jacobi reference frame for the coordinates and
velocities of the bodies, where the state vector of $m_{1}$ is defined as
$m_{0}$-centric, while the state vector of $m_{2}$ is measured with respect to
the barycenter of $m_{0}$ and $m_{1}$. We note that these coordinates and
velocities are required to compute the forcing terms given by
$\dot{\varphi_{i}}^{eq}$, $\epsilon_{i}^{\rho}$ and ${\epsilon_{i}}^{z}$, in
our differential equation system.

Regarding the resolution method, we choose the same path as the one taken in
Appendix 2 of Folonier et al. (2018). It consists in proposing a particular
stationary solution, inspired in the form (i.e. dependence on fundamental
frequencies) of the elliptical expansions of the forced terms
$\epsilon^{\rho}_{i}$ and $\dot{\varphi_{i}}^{eq}$. In this case, we propose a
solution of the form

$\displaystyle\Omega_{2}=\sum_{\@vec{\ell}}\\{\Omega_{2}\\}_{\@vec{\ell}}\cos{(l_{1}M_{1}+l_{2}M_{2}+l_{3}\varpi_{1}+l_{4}\varpi_{2}-\Phi_{\@vec{\ell},\Omega_{2}})}$
$\displaystyle
X_{2}=\sum_{\@vec{\ell}}\\{X_{2}\\}_{\@vec{\ell}}\cos{(l_{1}M_{1}+l_{2}M_{2}+l_{3}\varpi_{1}+l_{4}\varpi_{2}-\Phi_{\@vec{\ell},X_{2}})}$
$\displaystyle
Y_{2}=\sum_{\@vec{\ell}}\\{Y_{2}\\}_{\@vec{\ell}}\cos{(l_{1}M_{1}+l_{2}M_{2}+l_{3}\varpi_{1}+l_{4}\varpi_{2}-\Phi_{\@vec{\ell},Y_{2}})}$
(15)
$\displaystyle\mathcal{E}^{z}_{2}=\sum_{\@vec{\ell}}\\{\mathcal{E}^{z}_{2}\\}_{\@vec{\ell}}\cos{(l_{1}M_{1}+l_{2}M_{2}+l_{3}\varpi_{1}+l_{4}\varpi_{2}-\Phi_{\@vec{\ell},\mathcal{E}^{z}_{2}})},$

where $M_{1}$ and $\varpi_{1}$ are the mean anomaly and pericenter longitude
of the secondary star, while $M_{2}$ and $\varpi_{2}$ are those corresponding
to the planet. The amplitudes $\\{w\\}_{\@vec{\ell}}$ and the constant phases
$\Phi_{\@vec{\ell},x}$ depend on the subscripts
$\@vec{\ell}=(l_{1},l_{2},l_{3},l_{4})$ and each variable
$w=\Omega_{2},X_{2},Y_{2}$ or $\mathcal{E}^{z}_{2}$. This particular solution
and their derivatives are then introduced into the system (2.2) and (14), and
the coefficients can be calculated by equating terms with same trigonometric
argument and neglecting terms of higher-order. The terms in the system (4)
with $l_{1}=l_{2}=0$ correspond to the mean stationary solution and will be
given in what follows.

We choose to express the stationary solution of the rotational evolution
written in terms of the mean planetary spin $\langle\Omega_{2}\rangle$, mean
planetary lag-angle $\langle\delta_{2}\rangle$, mean planetary equatorial
flattening $\langle{\mathcal{E}}^{\rho}_{2}\rangle$ and mean planetary polar
flattening $\langle{\mathcal{E}}^{z}_{2}\rangle$. Up to 3rd-order in
$\alpha=a_{1}/a_{2}$ and to 2nd-order in the eccentricities $e_{1}$ and
$e_{2}$, including also the contribution of the eccentricities independent
4th-order and 6th-order terms, the solution is explicitly given by

$\displaystyle\langle\Omega_{2}\rangle$ $\displaystyle=$ $\displaystyle
n_{2}+6\frac{\gamma_{2}^{2}n_{2}}{\gamma_{2}^{2}+n_{2}^{2}}e_{2}^{2}+5\frac{\gamma_{2}^{2}n_{2}}{\gamma_{2}^{2}+n_{2}^{2}}e_{2}^{2}\mathcal{M}_{2}\alpha^{2}$
$\displaystyle-$
$\displaystyle\frac{325}{16}\frac{\gamma_{2}^{2}n_{2}}{\gamma_{2}^{2}+n_{2}^{2}}e_{1}e_{2}\cos(\Delta\varpi)\mathcal{M}_{3}\alpha^{3}$
$\displaystyle-$
$\displaystyle\frac{19}{2}\frac{\gamma_{2}^{2}(2n_{1}-2n_{2})}{\gamma_{2}^{2}+{(2n_{1}-2n_{2})}^{2}}\mathcal{M}_{2}^{2}\alpha^{4}+\langle\Omega_{2}\rangle_{o6}$
$\displaystyle\langle\delta_{2}\rangle$ $\displaystyle=$ $\displaystyle
6\frac{\gamma_{2}n_{2}}{\gamma_{2}^{2}+n_{2}^{2}}e_{2}^{2}+5\frac{\gamma_{2}n_{2}}{\gamma_{2}^{2}+n_{2}^{2}}e_{2}^{2}\mathcal{M}_{2}\alpha^{2}$
$\displaystyle-$
$\displaystyle\bigg{[}\frac{325}{16}\gamma_{2}\cos(\Delta\varpi)+\frac{175}{32}n_{2}\sin(\Delta\varpi)\bigg{]}\frac{n_{2}}{\gamma_{2}^{2}+n_{2}^{2}}e_{1}e_{2}\mathcal{M}_{3}\alpha^{3}$
$\displaystyle-$
$\displaystyle\frac{19}{2}\frac{\gamma_{2}(2n_{1}-2n_{2})}{\gamma_{2}^{2}+{(2n_{1}-2n_{2})}^{2}}\mathcal{M}_{2}^{2}\alpha^{4}+\langle\delta_{2}\rangle_{o6}$
$\displaystyle\frac{\langle{\mathcal{E}}^{\rho}_{2}\rangle}{\overline{\epsilon}^{\rho}_{2}}$
$\displaystyle=$ $\displaystyle
1+\bigg{(}\frac{3}{2}-4\frac{n_{2}^{2}}{\gamma_{2}^{2}+n_{2}^{2}}\bigg{)}e_{2}^{2}$
$\displaystyle+$
$\displaystyle\bigg{[}\frac{5}{4}+\frac{15}{8}e_{1}^{2}+\bigg{(}\frac{25}{4}-5\frac{n_{2}^{2}}{\gamma_{2}^{2}+n_{2}^{2}}\bigg{)}e_{2}^{2}\bigg{]}\mathcal{M}_{2}\alpha^{2}$
$\displaystyle-$
$\displaystyle\frac{25}{8}\frac{(3n_{2}^{2}+5\gamma_{2}^{2})}{\gamma_{2}^{2}+n_{2}^{2}}e_{1}e_{2}\cos(\Delta\varpi)\mathcal{M}_{3}\alpha^{3}$
$\displaystyle+$
$\displaystyle\bigg{[}\frac{105}{64}\mathcal{M}_{4}+\frac{4\gamma_{2}^{2}\mathcal{M}_{2}^{2}}{\gamma_{2}^{2}+{(2n_{1}-2n_{2})}^{2}}\bigg{]}\alpha^{4}+\frac{{\langle\mathcal{E}}^{\rho}_{2}\rangle_{o6}}{\overline{\epsilon}^{\rho}_{2}}$
$\displaystyle\frac{\langle{\mathcal{E}}^{z}_{2}\rangle}{\overline{\epsilon}^{M}_{2}}$
$\displaystyle=$
$\displaystyle\bigg{(}1+12\frac{\gamma_{2}^{2}}{\gamma_{2}^{2}+n_{2}^{2}}e_{2}^{2}+10\frac{\gamma_{2}^{2}}{\gamma_{2}^{2}+n_{2}^{2}}e_{2}^{2}\mathcal{M}_{2}\alpha^{2}$
$\displaystyle-$
$\displaystyle\frac{325}{8}\frac{\gamma_{2}^{2}}{\gamma_{2}^{2}+n_{2}^{2}}e_{1}e_{2}\cos(\Delta\varpi)\mathcal{M}_{3}\alpha^{3}$
$\displaystyle-$ $\displaystyle
19\frac{{(2n_{1}-2n_{2})}\gamma_{2}^{2}}{n_{2}(\gamma_{2}^{2}+{(2n_{1}-2n_{2})}^{2})}\mathcal{M}_{2}^{2}\alpha^{4}\bigg{)}$
$\displaystyle+$
$\displaystyle\frac{\overline{\epsilon}^{\rho}_{2}}{\overline{\epsilon}^{M}_{2}}\
\bigg{(}\frac{1}{2}+\frac{3}{4}e_{2}^{2}+\bigg{[}\frac{9}{8}+\frac{27}{16}e_{1}^{2}+\frac{45}{8}e_{2}^{2}\bigg{]}\mathcal{M}_{2}\alpha^{2}$
$\displaystyle-$
$\displaystyle\frac{225}{16}e_{1}e_{2}\cos(\Delta\varpi)\mathcal{M}_{3}\alpha^{3}+\frac{225}{128}\mathcal{M}_{4}\alpha^{4}\bigg{)}+\frac{{\langle\mathcal{E}}^{z}_{2}\rangle_{o6}}{\overline{\epsilon}^{M}_{2}}$

where $\langle\Omega_{2}\rangle_{o6}$, $\langle\delta_{2}\rangle_{o6}$,
${\langle\mathcal{E}}^{\rho}_{2}\rangle_{o6}$ and
${\langle\mathcal{E}}^{z}_{2}\rangle_{o6}$ are the terms independent on the
eccentricities, expanded up to 6th-order in $\alpha$ and listed in the
Appendix. We mention that the 5th-order terms are eccentricities dependent
and, thus, were excluded in this approximation. The mass factors are defined
in terms of the binary mass ratio $\beta=m_{1}/m_{0}$ as

$\displaystyle\mathcal{M}_{2}=\frac{\beta}{{(1+\beta)}^{2}}\ ,\
\mathcal{M}_{3}=\frac{\beta-\beta^{2}}{{(1+\beta)}^{3}}\ ,\
\mathcal{M}_{4}=\frac{\beta-\beta^{2}+\beta^{3}}{{(1+\beta)}^{4}}$ (20)

and the equilibrium mean flattenings are given by

$\displaystyle\overline{\epsilon}^{\rho}_{2}=\frac{15(m_{0}+m_{1})\mathcal{R}_{2}^{3}}{4m_{2}a^{3}_{2}}\
\ \ $ $\displaystyle\mathrm{and}$ $\displaystyle\ \ \
\overline{\epsilon}^{M}_{2}=\frac{5\mathcal{R}_{2}^{3}n_{2}^{2}}{4m_{2}\mathcal{G}}.$
(21)

Figure 4: Mean rotational speed $\langle\Omega_{2}\rangle$ (upper left panel),
mean lag-angle $\langle\delta_{2}\rangle$ (upper right panel), mean equatorial
flattening $\langle{\mathcal{E}}^{\rho}_{2}\rangle$ (lower left panel) and
mean polar flattening $\langle{\mathcal{E}}^{z}_{2}\rangle$ (lower right
panel) as a function of the planetary relaxation factor normalized by the mean
motion frequency $\gamma_{2}/n_{2}$. The initial conditions for the system
were taken from Table 1, but considering different planetary eccentricities
shown with different color curves, as indicated in the upper left panel.
Different types of curves represent different methods for obtaining the mean
values: full curves represent the numerical results, dashed curves represent
our 6th-order fit and dotted curves the results of considering the 2BP with a
central mass equal to the sum of the binary components. Vertical grey line in
the upper left panel represents the position in the diagram of a solid Earth
(Ferraz-Mello, 2013) and big color dots represent the final mean value
obtained with different $\gamma_{2}$-values, from simulations of Figure 3
captured in the 1:1 spin-orbit resonance: $\gamma_{2}/n_{2}=0.1$ (orange),
$\gamma_{2}/n_{2}=1$ (green), $\gamma_{2}/n_{2}=10$ (red) and
$\gamma_{2}/n_{2}=100$ (magenta).  Figure 5: Mean rotational speed
$\langle\Omega_{2}\rangle$ (first row panels), mean lag-angle
$\langle\delta_{2}\rangle$ (second row panels), mean equatorial flattening
$\langle{\mathcal{E}}^{\rho}_{2}\rangle$ (third row panels) and mean polar
flattening $\langle{\mathcal{E}}^{z}_{2}\rangle$ (last row panels) as a
function of the semimajor axes ratio $\alpha=a_{1}/a_{2}$, for the case
$\gamma_{2}=100n_{2}$. The initial conditions for the system were taken from
Table 1, but considering different planetary eccentricities shown with
different color curves (same color convention as in Figure 4) and different
secondary mass: $m_{1}=1/16m_{\odot}$ (first column) and $m_{1}=1/4m_{\odot}$
(second column). Different types of curves represent different methods for
obtaining the mean values with the same convention than in Figure 4 but
including the 4th-order fit with dotted-dashed curves. The magenta dots in the
right panels represent the final mean value obtained for this relaxation
factor from the simulations of Figure 3. Figure 6: Stationary spin map
predicted by our analytical model (equation (4)) for a CB planet in the
gaseous regime ($\gamma_{2}>>n_{2}$) and in an orbit with
$\alpha=a_{1}/a_{2}=1/3$, as a function of its eccentricity $e_{2}$ and the
binary mass ratio $m_{1}/m_{0}$. The white curve indicate the perfectly
synchronous state $\langle\Omega_{2}\rangle/n_{2}=1$. The white dots represent
the position in this diagram of the Kepler CB planets with expected large
gaseous component located close to the binary.

We first note from expressions (4), (4), (4) and (4) that the mean solutions
depend on the mean orbital configuration mainly through $e_{2}$ up to low-
order in $\alpha$, but also through $e_{1}$ and $\Delta\varpi$ up to higher-
orders. The dependence of the expressions with the masses $m_{0}$ and $m_{1}$,
and the semimajor axes ratio $\alpha$ are given in such a way that when
$m_{1}\to 0$ (i.e. a secondary star of very low mass) or $\alpha\to 0$ (i.e. a
planet far away from the binary), we recover the 2BP expressions of Folonier
et al. (2018).

On the other hand, according to preliminary numerical simulations, the mean
rotational solutions only exhibit a very weak dependence with $e_{1}$. This
fact is also expected if we notice that, in the mean analytical expression,
the functional dependence with this parameter only appears up to 3rd-order in
$\alpha$.

Another interesting feature of the mean solutions is the strong correlation
between the stationary spin $\langle\Omega_{2}\rangle$ with the lag-angle
$\langle\delta_{2}\rangle$: if we neglect the dependence on $\Delta\varpi$
(for example, for aligned and non-evolving pericenters or for the circulating
case), these quantities are simply related as $\langle\Omega_{2}\rangle-
n_{2}=\langle\delta_{2}\rangle/\gamma_{2}$. As a result, a sub or super-
synchronous mean stationary rotation speed is directly associated with a delay
or advance of the tidal bulge.

The rheology of the planet is taken into account through the ratio mainly
between the relaxation factor $\gamma_{2}$ and the main frequency of the
problem $n_{2}$. However, up to higher $\alpha$-order, another frequencies
such as $2n_{1}-2n_{2}$ starts become important in the solutions. For this
last reason, as we comment in Section 3, the exact position of the frontier
between the stiff regime and the gaseous regime is expected to be close to
$\gamma_{2}=n_{2}$, but not exactly at this point.

In Figure 4, we compare the mean values of $\langle\Omega_{2}\rangle$,
$\langle\delta_{2}\rangle$, $\langle{\mathcal{E}}^{\rho}_{2}\rangle$ and
$\langle{\mathcal{E}}^{z}_{2}\rangle$ obtained with different methods, as a
function of the normalized relaxation factor $\gamma_{2}/n_{2}$. The full
curves indicate the results of a numerical integration from an initial
condition far from the equilibrium which is allowed to evolve towards the
stationary state and then is filtered over the fast angles. The dashed curves
represent the results obtained respectively with equations (4) , (4), (4) and
(4), while the dotted curves represent the simplified 2BP case in which the
central binary is replaced by a unique body with mass equal to the sum of the
binary components. Big color dots represent the final mean value obtained with
different $\gamma_{2}$-values in simulations of Section 3, which result in a
capture in the 1:1 spin-orbit resonance (see Figure 3): the case of
$\gamma_{2}/n_{2}=0.1$ (in orange), $\gamma_{2}/n_{2}=1$ (in green),
$\gamma_{2}/n_{2}=10$ (in red) and $\gamma_{2}/n_{2}=100$ (in magenta). We
note that, in this last case, the magenta dots are not perfectly located on
the blue full curves due to this final condition, taken from simulation of
Figure 3, is almost but still not perfectly in the stationary state.

Figure 4 allows us to observe the existence of the two different regimes
commented in Section 3: the stiff regime for $\gamma_{2}<n_{2}$ and the
gaseous regime for $\gamma_{2}>n_{2}$, where far enough away from the
transition value $\sim n_{2}$, the rotational equilibrium solution becomes
independent of the planetary relaxation factor $\gamma_{2}$. In particular, in
the case of the mean stationary speed $\langle\Omega_{2}\rangle$ (upper left
panel), planets in the stiff regime are expected to asymptotically synchronize
their spins with the mean motion, independently of their orbital parameters,
body masses and viscosities. This prediction is completely analog to that
expected with the 2BP (Folonier et al., 2018). However, a much greater wealth
of behavior occurs in the gaseous regime: the stationary solutions do depend
on the orbital parameters and binary component masses. In this regime, unlike
what is expected in the 2BP (dotted curves), the presence of an inner
secondary star leads to lower stationary spins which can be sub-synchronous,
in the low planetary eccentric limit.

With respect to the accuracy of our analytical model in fitting the filtered
numerical integrations, we can observe that is not only able to predict the
sub-synchronous solution of the spins, but also gives a very precise
estimation of the equilibrium solution in the complete range of
$\gamma_{2}$-values. On the other hand, the 2BP simplification can also give a
good prediction, but the accuracy of the approximation is strongly dependent
on the initial condition and the relaxation factor range.

In Figure 5, we evaluate with more detail the accuracy of our mean stationary
solution given by expressions (4), (4), (4) and (4), as a function of the
small parameters $\alpha=a_{1}/a_{2}$ and $e_{2}$, for the gaseous case
$\gamma_{2}=100n_{2}$. Different $\alpha$-values were obtained by varying the
semi-major axis of the planet $a_{2}$. Although it is not a small parameter of
our model, we also explore a different secondary mass $m_{1}$.

We observe from Figure 5 that our 6th-order analytical fit is a very precise
estimation of the filtered numerical integrations for the rotational evolution
of a CB planet in a Kepler 38-like system. This occurs even for high
$\alpha$-values such as $\alpha\sim 0.4$, where the planet is probably inside
the instability limit imposed by the presence of an internal secondary
disturber (e.g. Holman & Wiegert, 1999). A good approximation is also obtained
with the 4th-order fit. However, because most Kepler CB planets are located
very close to binary (typically $\alpha\sim 0.3$), a higher order approach is
required to model more accurately its mean rotational behavior.

As expected, for the standard case ($m_{1}=1/4m_{\odot}$, right column) the
2BP simplification only is a good fit of the numerical simulation for low
$\alpha$-values, where the tidal effects are probably of little importance.
Curiously, the case of the mean polar flattening
$\langle{\mathcal{E}}^{z}_{2}\rangle$ presents an exception: the numerical
stationary solutions show little dependence with $\alpha$ and, while the 6th-
order fit is the best fit in the complete range, the 2BP simplification is
also a good approximation, even better than the 4th-order fit for planetary
positions very close to the binary. This may seem less strange if we remember
the expression of the equilibrium polar flattening $\epsilon_{i}^{z}$ given in
equation 2.1. Due to the definition of the flattenings chosen in this study
(see Folonier et al. (2018)), $\epsilon_{2}^{z}$ contains terms that depend on
the equatorial flattenings $\epsilon_{2,0}^{\rho}$ and $\epsilon_{2,1}^{\rho}$
which increase for higher $\alpha$-values, in addition to the classical
Maclaurin flattening $\epsilon_{2}^{M}$ which decreases as we move closer to
the binary, since the stationary planetary spin does. Thus, we expect a
similar behavior in the real polar flattening $\mathcal{E}_{2}^{z}$ that, due
to a potential balance of this two contributions, may result in a weak
dependence with $\alpha$.

The case of $m_{1}=1/16m_{\odot}$ (left panels of Figure 5) is considered as
an illustrative example that exhibits the strong dependence of the solution
with the secondary mass: for a low-mass secondary companion, the 2BP
simplification begins to be a good approximation, specially for CB planets in
wide orbits.

As an application example of our model predictions, in Figure 6 we show a
stationary spin map constructed with equation (4), for a CB planet in the
gaseous regime $\gamma_{2}/n_{2}=100$ located at $\alpha=1/3$, as a function
of its eccentricity $e_{2}$ and the binary mass ratio $m_{1}/m_{0}$. The white
curve represents the perfect synchronous state
$\langle\Omega_{2}\rangle=n_{2}$, which divides the map in a region where the
stationary spins are expected to be super-synchronous (above), and another one
where sub-synchronous spins are predicted (behind). With this reference, we
can observe that the size of the sub-synchronous region is proportional to the
normalized secondary mass, and tends to disappear for low secondary companion.
In competition with $m_{1}/m_{0}$, as we noticed in Zoppetti et al. (2019,
2020) and similar to that expected with the 2BP, the planetary eccentricity
$e_{2}$ leads to more rapid stationary spins.

We have included in the Figure 6 the nominal location in the diagram of the
Kepler CB planets observed very close to the binary and with an estimated mass
which allows us to infer that they have an important gaseous component
($m_{2}\gtrsim m_{Neptune}$): Kepler-16 (Doyle et al., 2011), Kepler-35 (Welsh
et al., 2012), Kepler-38 (Orosz et al., 2012), Kepler-47b (Orosz et al.,
2012), Kepler-64 (Schwamb et al., 2013), Kepler-413 (Kostov et al., 2014) and
Kepler-1661 (Socia et al., 2020). If we assume that all these CB planets have
an associated relaxation factor that satisfy $\gamma_{2}>>n_{2}$ and that, in
turn, these viscosities allowed them to have already reached their equilibrium
state, so the expected stationary spins for them is typically sub-synchronous.
The exceptional case is Kepler-413, an eccentric CB planet orbiting a binary
with two comparable-mass central stars, while Kepler-64 is expected to be very
close to the perfect synchronous state. However, as we mention in Section 3,
the amount of sub-synchronous shift is much less important than the expected
with the CTL tide model built in our previous works Zoppetti et al. (2019,
2020). The reasons of this difference will be explained in detail in the
following Section.

##  5 Revisiting the rotation of a CB planet in the CTL model

In the limit of gaseous bodies (i.e. $\gamma_{2}>>n_{1},n_{2}$), the
stationary spin equation (4) expanded up to 2nd-order in $\alpha$ becomes

$\displaystyle\langle\Omega_{2}\rangle_{st}^{gas}$ $\displaystyle=$
$\displaystyle
n_{2}+6n_{2}e_{2}^{2}+5n_{2}e_{2}^{2}\mathcal{M}_{2}\alpha^{2}.$ (22)

We notice that this expression is different than the one we obtained in
equation (22) of Zoppetti et al. (2019), using the classic CTL tidal model,
based on the expressions of the direct torques of Mignard (1979). This Section
is dedicated to investigate the origin of this discrepancy and comparing these
two different tidal models (i.e. CTL and creep) in the framework of the 3BP
and, thus, to obtain conclusions about the advantages and disadvantages of
each one, in different contexts.

In Section (2.1) of Zoppetti et al. (2019) we used a simplified model to test
the role of the cross tides in the 3BP with the CTL tidal model of Mignard
(1979). We considered only one extended mass $m_{0}$ under tidal interaction
of only one companion $m_{1}$ and used a third mass $m_{2}$ only as a test
particle, to evaluate the secular contribution of the cross torque, neglecting
the gravitational perturbations of $m_{1}$. We used the usual astrocentric
coordinate system, centered on the extended mass. In that simplified system we
found, employing numerical integrations and also analytical calculations, that
the secular net contribution of the cross torques on the test particle were
null. Then, we extrapolated that result to the general 3-extended body problem
in a Jacobian frame and neglected the cross forces and tides of Mignard
expressions. In what follows, we reformulate and improve the CTL model
developed in Zoppetti et al. (2019) with the inclusion of the cross tides, and
compare with our previous results.

Once again, we consider the 3BP composed of $m_{i}$, $m_{j}$ and $m_{k}$
represented in Figure 7, where between each pair of bodies we compute the
complete tidal interaction, in addition to the gravitational forces. As a
result of this consideration of the tides, due to the presence of each body,
let’s say $m_{i}$, one ellipsoid is generated on $m_{j}$ and one ellipsoid on
$m_{k}$, both represented as grey ellipsoids on each respective mass of Figure
7.

Figure 7: Sketch of the 3BP considering all the tidal interactions between
pairs of bodies. The grey ellipsoids are induced on the companions due to the
presence of $m_{i}$, the red ellipsoids due to the presence of $m_{j}$ and the
blue ellipsoids are due to the presence of $m_{k}$.

Let us denote as ${\bf{F}}_{i,j}^{k}$ to the tidal force applied on $m_{i}$
due to the ellipsoid that the presence of $m_{k}$ generated on $m_{j}$. In
Figure 7, ${\bf{F}}_{i,j}^{k}$ corresponds to the force applied on $m_{i}$ due
to the blue ellipsoid generated on $m_{j}$. Thus, if we take into account all
the action forces applied on $m_{i}$, we must consider the tides caused by the
two ellipsoids that $m_{i}$ itself generates on its companions (grey
ellipsoids) plus the tides due to the cross ellipsoids (blue ellipsoid on
$m_{j}$ and red ellipsoid on $m_{k}$). On the other hand, to complete the
force diagram, we also need to consider the reaction forces that $m_{i}$
exerts on its two companions, due to the ellipsoids generated on itself (blue
and red ellipsoids on $m_{i}$). As a result, the total tidal force applied on
$m_{i}$ is given by

$\displaystyle{\bf{F}}_{i}=\underbrace{({\bf{F}}_{i,j}^{i}+{\bf{F}}_{i,j}^{k}+{\bf{F}}_{i,k}^{i}+{\bf{F}}_{i,k}^{j})}_{action\
tides}-\underbrace{({\bf{F}}_{j,i}^{j}+{\bf{F}}_{j,i}^{k}+{\bf{F}}_{k,i}^{k}+{\bf{F}}_{k,i}^{j})}_{reaction\
tides}$ (23)

We note that the forces with three different indexes in equation (23)
correspond to the cross tides, which were neglected in Zoppetti et al. (2019),
while the forces with two repeated indexes correspond to the direct tides,
which we did consider in our previous work.

If we adopt the classical formalism of Mignard (1979) to model the tidal
interactions between each pair of bodies with, the explicit expression of the
force ${\bf{F}}_{i,j}^{i}$ is expanded up to first order in the time lag
$\Delta t_{j}$ as

${\bf{F}}_{i,j}^{i}={\bf{F}}_{i,j}^{k,(0)}+\Delta
t_{j}{\bf{F}}_{i,j}^{k,(1)},$ (24)

where (see Zoppetti et al. (2019))

$\displaystyle{\bf{F}}_{i,j}^{k,(0)}=\frac{3}{2}\mathcal{G}\mathcal{R}_{j}^{5}k_{2,j}\frac{m_{i}}{\Delta_{ij}^{5}}\frac{m_{k}}{\Delta_{kj}^{5}}\bigg{[}2({\bf\Delta}_{ij}\cdot{\bf\Delta}_{kj}){\bf\Delta}_{kj}+$
$\displaystyle(\Delta_{kj}^{2}-\frac{5}{\Delta_{ij}^{2}}{({\bf\Delta}_{ij}\cdot{\bf\Delta}_{kj})}^{2}){\bf\Delta}_{ij}\bigg{]},$
(25)
$\displaystyle{\bf{F}}_{i,j}^{k,(1)}=3\mathcal{G}\mathcal{R}_{j}^{5}k_{2,j}\frac{m_{i}}{\Delta_{ij}^{5}}\frac{m_{k}}{\Delta_{kj}^{5}}\bigg{[}\frac{({\bf\Delta}_{kj}\cdot\dot{{\bf\Delta}_{kj}})}{\Delta_{kj}^{2}}\bigg{(}5{\bf\Delta}_{kj}({\bf\Delta}_{ij}\cdot{\bf\Delta}_{kj})-\Delta_{kj}^{2}{\bf\Delta}_{ij}\bigg{)}$
$\displaystyle-[{\bf\Delta}_{kj}\cdot({\bf\Omega}_{j}\times{\bf\Delta}_{ij})+{\bf\Delta}_{ij}\cdot\dot{{\bf\Delta}_{kj}}]{\bf\Delta}_{kj}-[{\bf\Delta}_{kj}\times{\bf\Omega}_{j}+\dot{{\bf\Delta}_{kj}}]({\bf\Delta}_{ij}\cdot{\bf\Delta}_{kj})$
$\displaystyle+\frac{5}{\Delta_{ij}^{2}}{\bf\Delta}_{ij}\bigg{(}({\bf\Delta}_{ij}\cdot{\bf\Delta}_{kj})[{\bf\Delta}_{kj}\cdot({\bf\Omega}_{j}\times{\bf\Delta}_{ij})+{\bf\Delta}_{ij}\cdot\dot{{\bf\Delta}_{kj}}]$
$\displaystyle-\frac{({\bf\Delta}_{kj}\cdot\dot{{\bf\Delta}_{kj}})}{2\Delta_{kj}^{2}}(5{({\bf\Delta}_{ij}\cdot{\bf\Delta}_{kj})}^{2}-\Delta_{ij}^{2}\Delta_{kj}^{2})\bigg{)}\bigg{]},$
(26)

with $k_{2,j}$ the second degree Love number of $m_{j}$.

The conservation of the total angular momentum of the system imposes that the
total orbital angular momentum ${\bf{L}}_{orb}$ must satisfy

$\sum_{i=0}^{2}C_{i}\dot{\bf{\Omega}_{i}}=-\dot{\bf{L}}_{orb},$ (27)

where the variation of the total orbital momentum can be written in terms of
the relative positions as (Zoppetti et al., 2019)

$\displaystyle\dot{\bf{L}}_{orb}={\bf\Delta}_{j,i}\times({\bf{F}}_{j,i}^{j}+{\bf{F}}_{j,i}^{k})+{\bf\Delta}_{k,i}\times({\bf{F}}_{k,i}^{k}+{\bf{F}}_{k,i}^{j})+$
(28)
$\displaystyle{\bf\Delta}_{i,j}\times({\bf{F}}_{i,j}^{i}+{\bf{F}}_{i,j}^{k})+{\bf\Delta}_{k,j}\times({\bf{F}}_{k,j}^{k}+{\bf{F}}_{k,j}^{i})+$
$\displaystyle{\bf\Delta}_{i,k}\times({\bf{F}}_{i,k}^{i}+{\bf{F}}_{i,k}^{j})+{\bf\Delta}_{j,k}\times({\bf{F}}_{j,k}^{j}+{\bf{F}}_{j,k}^{i})$

As discussed in Section 4, equation (27) may be decoupled if we assume that
the variation in the spin angular momentum of the body $m_{j}$ is only due to
the terms in the equation associated with its own deformation (see Zoppetti et
al. (2019)). Using this boundary condition, we have

$\displaystyle
C_{i}\dot{\bf{\Omega}}_{i}=-\bigg{[}\underbrace{({\bf\Delta}_{ji}\times{\bf{F}}_{ji}^{j}+{\bf\Delta}_{ki}\times{\bf{F}}_{ki}^{k})}_{direct\
torques}+\underbrace{({\bf\Delta}_{ji}\times{\bf{F}}_{ji}^{k}+{\bf\Delta}_{ki}\times{\bf{F}}_{ki}^{j})}_{cross\
torques}\bigg{]}$ (29)

Again, we mention that the rotational evolution descripted by equation (15) of
Zoppetti et al. (2019), only considers the direct torques, i.e. the first two
terms in equation (29).

Figure 8:  Rotational evolution of a CB planet with the CTL model of Mignard
(1979) (equation (29)). The initial conditions were taken from Table 1, with
the exception of the initial planetary eccentricity ($e_{2}(0)=0.05$ in the
upper panel and $e_{2}(0)=0.15$ in the bottom panel). Full curves indicate the
results of filtered numerical integration: excluding the cross tides (in
yellow) and considering the full expression (in magenta). Dashed curves
represents the results of different analytical models: our previous CTL tidal
model (Zoppetti et al., 2019) (in brown) and the results of the analytical
creep tidal model presented in equation (4) of this work, expanded up to
2-order in $\alpha$ and considering the gaseous limit
($\gamma_{2}>>n_{1},n_{2}$), are shown in grey.

In Figure 8, we compare the different rotational evolution that results from
the CTL tidal model of Mignard (1979) in the 3BP, by considering the full
expression in equation (29) and the version constructed neglecting the cross
tides. The initial conditions were taken from Table 1, except for the initial
planetary eccentricity which was varied: $e_{2}(0)=0.05$ in the upper panel
and $e_{2}(0)=0.15$ in the bottom panel. In full curves we show the results of
digitally filtered numerical integrations of equation (29), considering the
full expression of the torque equation (in magenta) and neglecting the cross
tides (in yellow). We have also included in this simulations the punctual
gravitational interactions between the bodies, to compare with our previous
simulations of Zoppetti et al. (2019).

The results obtained with different analytical models are shown with dashed
curves in Figure 8: with brown curves we represent the results of the
analytical CTL tidal model of Zoppetti et al. (2019) (equation (22)) and with
grey curves we represent the results of the analytical creep tidal model
presented in this work (equation (4)) expanded up to 2-order in $\alpha$, in
the gaseous case ($\gamma_{2}>>n_{1},n_{2}$).

The comparison between the numerical integrations of the torque equation (29)
without cross tides and the analytical results of the CTL model of Zoppetti et
al. (2019) are analogous to the one presented in the right panels of Figure 7
of that work. However, here we adopt a tidal parameter
$Q_{2}^{\prime}=3/(2n_{2}k_{2,2}\Delta t_{2})$ that has been chosen in order
to the associated relaxation factor satisfy $\gamma_{2}=100n_{2}$ (see Ferraz-
Mello (2013)) and, thus, the rotational evolution timescales be similar to the
one presented in the magenta simulation of Figure 3.

We first notice that the stationary solution of $\langle\Omega_{2}\rangle$
exhibits long-term periodic oscillations, specially for the initially
eccentric case. As mention in Zoppetti et al. (2019), this is due to the
secular variation induced on the planetary eccentricity by the gravitational
punctual interactions.

Regarding the comparison with the numerical simulations, we can notice in both
panels of Figure 8, that the analytical model presented in Zoppetti et al.
(2019) (brown curves) fits very well the stationary value of the planetary
rotation, when we neglect the cross tides in the torque equation. This is an
expected result since, in the construction of that model, we took into account
only the direct tidal forces and not the general expression of equation (5).

The comparison of the numerical integration for the full torque equation with
the analytical results obtained with the creep model (in the gaseous limit),
reveals two very important facts. On one hand, from the proximity of the
magenta curves with the grey curves, we show that the results obtained for the
rotational stationary state with the CTL tidal model of Mignard (1979) are the
same that the one obtained with the creep model in the limit of very large
relaxation factor ($\gamma_{2}>>n_{1},n_{2}$). This result was checked in a
completely analytical calculation, by expanding the position and velocity
vectors of equation (29) in $\alpha$ and the eccentricities, and then
averaging with respect to both mean longitudes, with the help of an algebraic
manipulator.

On the other hand, the difference between the magenta and yellow curves of
Figure 8, shows that, unlike what we assumed in Zoppetti et al. (2019) from
analyzing the simple case in which we considered only one extended mass and a
reference frame in solidarity with it, the cross tides are not negligible in
the general self-consistent problem of three extended masses. However, the
physical effect of the tides due to the presence that a binary (instead a
single star) exerts on an external planet is the same as we reported in our
previous works: it shifts the rotational stationary solution to lower spin
values, for an amount that is proportional to the semimajor axes ratio
$\alpha$ and the secondary mass $m_{1}$. The no-cross tides model developed in
Zoppetti et al. (2019, 2020) simply overestimate this effect.

##  6 Summary and discussion

In this work, we present a model for treating the tides in the general 3BP, in
the case in which all the bodies are considered extended and tidally
interacting. In particular, we consider the planar case where all the bodies
have spins vectors pointing in a direction perpendicular to the orbital plane.
We use the creep tide formalism first presented in Ferraz-Mello (2013) to
compute the tidal forces and torques, and apply it to study the rotational
evolution of a circumbinary planet. The polar flattening due to the rotation
of the bodies is also naturally taken into account in the theory (Folonier et
al., 2018).

First, we present the creep differential equations that govern the rotational
evolution of the extended bodies, in the framework of the 3BP. We show that,
on each body, the tidal interactions due to the presence of the two companions
can be modelled by an equivalent unique interaction. This allows the formalism
to be extended to the general N-bodies tidally interacting case, where the
difficulty lies in finding the parameters of the equivalent interaction. As in
the 2BP, the physical quantities that characterize the rotational evolution of
each body in this model are the spin velocities, the equivalent lag-angles,
tidal flattenings and polar flattenings.

We then apply the model to study the rotational evolution of a circumbinary
planet, using as a working example a Kepler-38 like system (Orosz et al.,
2012). To do this, we take advantage of the adiabatic nature of the rotational
respect to the orbital evolution in this problem, which let us assume as fixed
the orbital elements with the exception of the fast angles.

Using direct numerical integrations, we find that, for planets in low
eccentric orbits and viscosities in the range of the ones estimated for Solar
System planets by Ferraz-Mello (2013), the more probably final stationary
state is the 1:1 spin-orbit resonance, with the equivalent lag-angle
oscillating around 0. The required timescales for reaching the stationary
state depends on the relaxation factor of the bodies: taking as a reference
the relaxation factors estimated for Ferraz-Mello (2013), we expect that CB
stiff planets in a Kepler-38 like system synchronize its spins in timescales
much lower than the age of the system, while this is not necessary the case
for planets with a large gaseous component. However, the great uncertainty
about the internal structure of exoplanets and, even more, about their
associated viscosities inhibit an accurate estimate of its rotational
evolution timescale.

The high probability of capture in the 1:1 spin-orbit resonance, for low
eccentric CB planets with arbitrary viscosities, motivated us to study this
particular solution in more detail. With this aim, we search for a high-order
analytical expression for the mean solution of the rotational stationary
state. The extreme proximity to their central binary to which Kepler’s CB
planets have been observed, in addition to the strong dependence of tidal
effects with the relative distance between the bodies, forced us to develop a
high-order model in semimajor axes ratio. As a result of this, our solution
has shown to be very precise in predicting the full numerical integrations,
for arbitrary planetary viscosity and very high $\alpha$-values, even inside
the expected gravitationally instable region.

From a practical point of view, when studying long-term tidal evolution,
having a secular approximation of the stationary rotational solution will
considerably accelerate the numerical simulation of the full system, avoiding
having to integrate the rotational evolution and, instead, introducing an
analytical recipe for the equilibrium state. This semi-analytic method is
particularly useful in the case of working with gaseous bodies with large
relaxation factors, for which the creep equation becomes particularly complex
to solve efficiently.

As an application of our analytical model, we consider the Kepler CB planets
observed close to the central binary and with an expected large gaseous
component. If we assume that all these CB planets have associated viscosities
that let them to have reached their rotational steady state, our model
predicts that their typical rotation speeds are sub-synchronous. This result
is mainly due to the presence of a central binary (instead of a single star)
and is predicted by our model for very low eccentric CB planets close to a
binary with a high mass ratio. However, the stationary rotation for gaseous CB
planets predicted by the creep tide theory presented here is different than
the one expected with the CTL model presented in Zoppetti et al. (2019, 2020).

To explain this discrepancy, we finally discuss the results obtained for the
stationary spin of a CB planet, with the creep tide theory (in the gaseous
limit) and compare with those obtained using the CTL tidal model in Zoppetti
et al. (2019, 2020). We show that both model predict exactly the same solution
as long as the cross torques are taken into account in the CTL formalism.
These torques had been neglected in our previous works, from the study of the
classical case of a single extended body with a system of coordinates attached
to it. However, for the general case of extended 3BP and an arbitrary
reference frame, this is not necessarily true. Actually, one of the main
advantage of the creep tide theory over CTL model is the treatment of cross
torques. While in the CTL model the general expressions for the tidal force
and torques are very complex and only valid for small and constant time-lags,
in the creep theory these quantities are calculated from the equilibrium
figures where the cross torques are naturally taken into account.

In addition to this, the CTL model seems to be contained within the creep tide
theory, for the case of relaxation factors much higher than the fundamental
frequencies of the system. Moreover, in the creep tide formalism the lag-angle
is not an external ad-hoc quantity which needs to be small and constant, but
can take large values and its temporal evolution is solved from a differential
equation. Finally, among other of its most interesting advantages, we name
that the calculation of an equivalent equilibrium figure allows not only to
take into account the cross torques in a natural way but also to incorporate
some additional physical phenomena such as polar flattening. As a result, the
theory also provide the evolution of the geometric shape of the interacting
bodies.

Among its main limitations, we found some difficulties in numerically solving
the creep equation for gaseous bodies with high relaxation factors in an
efficient way. This means that, when studying long-term tidal evolution, some
analytical recipes need to be incorporated to accelerate the numerical
integrations, such as the mean values that describe the rotational evolution
presented in this work. On the other hand, as noticed in Ferraz-Mello (2015),
the creep theory is analog to those which consider as Maxwell bodies to the
extended masses (Correia et al., 2014), but without the elastic component of
the tide. More complex rheologies such as Sundberg-Cooper model (Sundberg &
Cooper, 2010) may be neccesary to reproduce many of the intricacies seen in
the deformation of real materials (e.g. Renaud & Henning, 2018).

## Appendix: 6th-order terms for the mean stationary solution

The non-eccentric terms expanded up to 6th-order in $\alpha$ of expressions
(4), (4), (4) and (4) are respectively given by

$\displaystyle\langle\Omega_{2}\rangle_{o6}=$
$\displaystyle-\frac{5\mathcal{M}_{2}(n_{1}-n_{2})\gamma_{2}^{2}\alpha^{6}}{32\left(\gamma_{2}^{2}+(n_{1}-n_{2})^{2}\right)\left(\gamma_{2}^{2}+(2n_{1}-2n_{2})^{2}\right)\left(\gamma_{2}^{2}+(3n_{1}-3n_{2})^{2}\right)}$
(30) $\displaystyle\times$ $\displaystyle\bigg{[}5\
\bigg{(}121\gamma_{2}^{4}+805\gamma_{2}^{2}(n_{1}-n_{2})^{2}+744(n_{1}-n_{2})^{4}\bigg{)}\mathcal{M}_{4}$
$\displaystyle-\bigg{(}729\gamma_{2}^{4}+5265\gamma_{2}^{2}(n_{1}-n_{2})^{2}+4836(n_{1}-n_{2})^{4}\bigg{)}\mathcal{M}_{2}^{2}\bigg{]}$
$\displaystyle\langle\delta_{2}\rangle_{o6}=\frac{\langle\Omega_{2}\rangle_{o6}}{\gamma_{2}}$
(31)
$\displaystyle\frac{\langle\mathcal{E}^{\rho}_{2}\rangle_{o6}}{\overline{\epsilon}^{\rho}_{2}}=$
$\displaystyle\frac{5M_{2}\alpha^{6}}{256\left(\gamma_{2}^{2}+(n_{1}-n_{2})^{2}\right)\left(\gamma_{2}^{2}+4(n_{1}-n_{2})^{2}\right)\left(\gamma_{2}^{2}+9(n_{1}-n_{2})^{2}\right)}$
(32) $\displaystyle\times$
$\displaystyle\Big{[}(105\mathcal{F}_{1}+776\mathcal{F}_{2})\gamma_{2}^{6}+70(21\mathcal{F}_{1}+76\mathcal{F}_{2})\gamma_{2}^{4}(n_{1}-n_{2})^{2}$
$\displaystyle+(5145\mathcal{F}_{1}+5024\mathcal{F}_{2})\gamma_{2}^{2}(n_{1}-n_{2})^{4}+378\mathcal{F}_{1}(n_{1}-n_{2})^{6}\Big{]}$
$\displaystyle\frac{\langle{\mathcal{E}}^{z}_{2}\rangle_{o6}}{\overline{\epsilon}^{M}_{2}}=\frac{2\langle\Omega_{2}\rangle_{o6}}{n_{2}}+\frac{\overline{\epsilon}^{\rho}_{2}}{\overline{\epsilon}^{M}_{2}}\frac{1225}{512}\mathcal{M}_{2}\mathcal{F}_{1}\alpha^{6}$
(33)

where

$\displaystyle\mathcal{F}_{1}=\frac{1+\beta^{5}}{(1+\beta)^{5}}\ \ $ ;
$\displaystyle\ \ \mathcal{F}_{2}=\frac{\beta(1-\beta)^{2}}{(1+\beta)^{4}}.$
(34)

###### Acknowledgements.

This research was funded by CONICET, SECYT/UNC, FONCYT and FAPESP (Grant
2016/20189-9 and 2016/13750-6. The authors also wish to thank the anonymous
referee for helpful suggestions that greatly improved this work)

## References

  * Correia et al. (2014) Correia, A. C. M., Boué, G., Laskar, J., et al. 2014, A&A, 571, A50.
  * Correia et al. (2016) Correia, A. C. M., Boué, G., & Laskar, J. 2016, Celestial Mechanics and Dynamical Astronomy, 126, 189.
  * Doyle et al. (2011) Doyle, L. R., Carter, J. A., Fabrycky, D. C., et al. 2011, Science, 333, 1602.
  * Efroimsky (2012) Efroimsky, M. 2012, ApJ, 746, 150.
  * Efroimsky (2015) Efroimsky, M. 2015, AJ, 150, 98.
  * Ferraz-Mello et al. (2008) Ferraz-Mello, S., Rodríguez, A., & Hussmann, H. 2008, Celestial Mechanics and Dynamical Astronomy, 101, 171.
  * Ferraz-Mello (2013) Ferraz-Mello, S. 2013, Celestial Mechanics and Dynamical Astronomy, 116, 109.
  * Ferraz-Mello (2015) Ferraz-Mello, S. 2015, Celestial Mechanics and Dynamical Astronomy, 122, 359.
  * Ferraz-Mello (2015) Ferraz-Mello, S. 2015, A&A, 579, A97.
  * Folonier et al. (2015) Folonier, H. A., Ferraz-Mello, S., & Kholshevnikov, K. V. 2015, Celestial Mechanics and Dynamical Astronomy, 122, 183.
  * Folonier et al. (2018) Folonier, H. A., Ferraz-Mello, S., & Andrade-Ines, E. 2018, Celestial Mechanics and Dynamical Astronomy, 130, 78.
  * Gomes et al. (2019) Gomes, G. O., Folonier, H. A., & Ferraz-Mello, S. 2019, Celestial Mechanics and Dynamical Astronomy, 131, 56.
  * Holman & Wiegert (1999) Holman, M. J. & Wiegert, P. A. 1999, AJ, 117, 621.
  * Hurford et al. (2007) Hurford, T. A., Sarid, A. R., & Greenberg, R. 2007, Icarus, 186, 218.
  * Hut (1981) Hut, P. 1981, A&A, 99, 126.
  * Kostov et al. (2014) Kostov, V. B., McCullough, P. R., Carter, J. A., et al. 2014, ApJ, 784, 14.
  * Lambert (1991) Lambert, J. D. 1991, Numerical methods for ordinary differential systems: the initial value problem. John Wiley & Sons, Inc.
  * Makarov & Efroimsky (2013) Makarov, V. V. & Efroimsky, M. 2013, ApJ, 764, 27.
  * Mignard (1979) Mignard, F. 1979, Moon and Planets, 20, 301.
  * Mills & Mazeh (2017) Mills, S. M. & Mazeh, T. 2017, ApJ, 839, L8.
  * Murray & Dermott (1999) Murray, C. D. & Dermott, S. F. 1999, Solar system dynamics by C.D. Murray and S.F. McDermott. (Cambridge, UK: Cambridge University Press), ISBN 0-521-57295-9 (hc.), ISBN 0-521-57297-4 (pbk.).
  * Orosz et al. (2012) Orosz, J. A., Welsh, W. F., Carter, J. A., et al. 2012, ApJ, 758, 87.
  * Orosz et al. (2012) Orosz, J. A., Welsh, W. F., Carter, J. A., et al. 2012, Science, 337, 1511.
  * Renaud & Henning (2018) Renaud, J. P. & Henning, W. G. 2018, ApJ, 857, 98.
  * Socia et al. (2020) Socia, Q. J., Welsh, W. F., Orosz, J. A., et al. 2020, AJ, 159, 94.
  * Schwamb et al. (2013) Schwamb, M. E., Orosz, J. A., Carter, J. A., et al. 2013, ApJ, 768, 127.
  * Sundberg & Cooper (2010) Sundberg, M. & Cooper, R. F. 2010, Philosophical Magazine, 90, 2817.
  * Tisserand (1891) Tisserand, Félix. Traité de mécanique céleste. Gauthier-Villars, 1891.
  * Welsh et al. (2012) Welsh, W. F., Orosz, J. A., Carter, J. A., et al. 2012, Nature, 481, 475.
  * Zoppetti et al. (2018) Zoppetti, F. A., Beaugé, C., & Leiva, A. M. 2018, MNRAS, 477, 5301.
  * Zoppetti et al. (2019) Zoppetti, F. A., Beaugé, C., Leiva, A. M., et al. 2019, A&A, 627, A109.
  * Zoppetti et al. (2020) Zoppetti, F. A., Leiva, A. M., & Beaugé, C. 2020, A&A, 634, A12.