<table><caption>Table 5.1 — Value of $\kappa_2/Q$ for eight bodies in the Solar System. For the planets and satellites the values are constrained by astrometric measurements (Lainey, 2016). For the Sun, we used $\kappa_2 = 0.02$ (Claret and Cunha, 1997, Fig. 3) and $Q = 2 \times 10^5$ (Barker, 2020, Fig. 8).</caption><thead><tr><th></th><th>Earth</th><th>Moon</th><th>Mars</th><th>Jupiter</th><th>Io</th><th>Saturn</th><th>Uranus</th><th>Sun</th></tr></thead><tbody><tr><th>10<sup>4</sup> &kappa;<sub>2</sub>/Q</th><td>11</td><td>6.4</td><td>17</td><td>0.11</td><td>150</td><td>1.6</td><td>&gt; 2</td><td>0.001</td></tr></tbody></table>

## 5.3 The pseudo-Hamiltonian formalism

Combining together Eqs. (5.4), (5.12), (5.13) and (5.36), the final form for the perturbed potential due to the tidal bulge raised upon $\mathbb{B}$ is

$$V(\mathbf{r}) = -\kappa_2 \frac{G \tilde{m}}{r} \left(\frac{R_0}{r}\right)^2 \left(\frac{R_0}{r^*}\right)^3 P_2(\cos S^*), \quad (5.45)$$

where $\mathbf{r}^* = \tilde{\mathbf{r}}(t - \Delta t)$ and $S^*$ is the angle between $\mathbf{r}$ and $\mathbf{r}^*$. Studying a dynamical system with tidal dissipation, we write the Hamiltonian of the problem as

$$\mathcal{H} = \mathcal{H}_K(\mathbf{p}) + \iota_1 \mathcal{H}_P(\mathbf{p}, \mathbf{q}) + \iota_2 \mathcal{U}(\mathbf{p}, \mathbf{q}, \mathbf{p}^*, \mathbf{q}^*) + \mathcal{T}(\mathbf{p}), \quad (5.46)$$

where the tidal potential $\mathcal{U}$ contains every tidal contribution and $\mathcal{T}$ is the sum of the kinetic energies of rotation of the non point-mass bodies. The Hamiltonians $\mathcal{H}_K$ and $\mathcal{H}_P$ are the Keplerian and perturbative contributions of the point-mass Hamiltonian, while the vector $\mathbf{p}$ contains the momenta conjugated to the generalised coordinates $\mathbf{q}$. The order of magnitude of the parameter $\iota_1$ is the planet-to-star mass ratio, or the satellite-to-planet mass ratio, depending on the system, while $\iota_2$ is generally much smaller than $\iota_1$, due to the factor $R_0^5/r^5$ in Eq. (5.45). In the $N$-body problem, the tidal potential $\mathcal{U}$ is the sum of $N(N-1)^2$ contributions<fnref n="4" /> like in Eq. (5.45), but suitable simplifications allow most of them to be discarded, as we show in Sect. 6.1.1.

The equations of motion are obtained from $\mathcal{H}$ using the Hamilton Eqs. (2.18), considering that the quantities $\mathbf{p}^*$ and $\mathbf{q}^*$ are constant parameters; that is, we do not derive with respect to the starred variables. The subsequent vector field then depends on the parameters $\mathbf{p}^*$ and $\mathbf{q}^*$, and the dependency is removed with the first-order Taylor expansion

$$\mathbf{p}^* = \mathbf{p} - \Delta t \dot{\mathbf{p}}, \quad \mathbf{q}^* = \mathbf{q} - \Delta t \dot{\mathbf{q}}, \quad (5.47)$$

in the vector field. The Hamilton equations are used again in Eq. (5.47) to compute $\dot{\mathbf{p}}$ and $\dot{\mathbf{q}}$. To prevent an implicit relation, only the main contributions of the Hamiltonian are considered. In other words, we assume $\iota_1 = \iota_2 = 0$ to compute the time-derivatives in Eq. (5.47), that is

$$\mathbf{p}^* = \mathbf{p} \quad \text{and} \quad \mathbf{q}^* = \mathbf{q} - \Delta t \frac{\partial (\mathcal{H}_K + \mathcal{T})}{\partial \mathbf{p}}. \quad (5.48)$$

<footnote n="4">⁴Each of the $N$ bodies raises $N-1$ bulges, and every bulge affects the orbit of $N-1$ bodies.</footnote>
# Chapter 6

## Tides in the planar three-body co-orbital problem

*Most results of this chapter were first published in Couturier et al. (2021), except those of Sect. 6.4, which are original.*

### 6.1 Pseudo-Hamiltonian and equations of motion

This chapter extends Chap. 3 by adding tidal dissipation to the analytical model. We will add terms to the conservative Hamiltonian $\mathcal{H}_K + \mathcal{H}^{(0)} + \mathcal{H}^{(2)} + \mathcal{H}^{(4)}$ (Eqs. (3.24) & (3.29)) to account for tidal effects, using the results from Chap. 5. In this chapter, like in Chap. 3, $\iota$ denotes $(m_1 + m_2)/m_0$. The central body and two co-orbitals have radii $R_0$, $R_1$ and $R_2$, respectively. We showed in Sect. 5.2.3 that in order to have simple analytical expressions for the tidal perturbations, we need to consider the constant-$\Delta t$ model. We adopt this model in this chapter and the rheology of the bodies is hence characterized by their second Love number $\kappa_2$ (Eqs. (5.35) & (5.36)) and their quality factor $Q$ (Sect. 5.2.4).

#### 6.1.1 Expression of the pseudo-Hamiltonian

We want to build a tidal pseudo-Hamiltonian by adding to the conservative Hamiltonian contributions of the form (5.45). However, each of the three bodies can be a perturbing body, a perturbed body, or an interacting one, and taking into account every single tidal contribution will generate 12 new terms like in Eq. (5.45) (Sect. 5.3). In order to prevent a gruesome expression for the tidal pseudo-Hamiltonian, we show that for a star–planet–planet system, 10 out of 12 of these contributions can be neglected.

Indeed, tides raised by a planet on a planet, as well as tides raised by the star on a planet, and interacted with by the other planet, can readily be discarded due to the low masses and radii of the planets, with respect to that of the star. Hence, only six contributions remain. Two of these correspond to tides raised by the star on each planet
and interacted with by the star, and four correspond to tides raised by a planet on the star and interacted with by each planet. The four latter have an order of magnitude $\kappa_2^{(0)}/Q_0(\iota m_0)^2(R_0/\bar{a})^5$, while the two former have an order of magnitude $\kappa_2^{(j)}/Q_j m_0^2(R_j/\bar{a})^5$, where $j \in \{1, 2\}$. Assuming a similar density for all the bodies, we have $R_j \propto m_j^{1/3}$, and then (Table 5.1)

$$ \frac{\kappa_2^{(0)}(\iota m_0)^2 R_0^5 / Q_0}{\kappa_2^{(j)} m_0^2 R_j^5 / Q_j} = \iota^{1/3} \frac{\kappa_2^{(0)}}{\kappa_2^{(j)}} / Q_0 \ll 1. \qquad (6.1) $$

Tides raised by the planets on the star can be discarded too and we only consider the two contributions due to tides raised by the star on each planet and interacted with by the star. In the case of satellites co-orbiting a planet, we stress that tides raised on the central body cannot be neglected, because the $\kappa_2/Q$ of planet is not smaller than that of satellites, and the ratio (6.1) could be close to unity. We will treat this case in Sect. 6.4. For now, we write the perturbation to the conservative Hamiltonian, due to tides, as

$$ H_t = U_t^{(1)} + U_t^{(2)} + T_1 + T_2, \qquad (6.2) $$

where $T_j$ is the kinetic energy of rotation of planet $j$ and $U_t^{(j)}$ is the contribution due to tides raised by the star on planet $j$ and interacted with by the star. In the heliocentric reference frame, we have (Eq. (5.45))

$$ U_t^{(j)} = -\kappa_2^{(j)} \mathcal{G} m_0^2 \frac{R_j^5}{r_j^3 r_j^{*3}} P_2(\cos S) \quad \text{and} \quad T_j = \frac{\Theta_j^{\prime 2}}{2C_j}, \qquad (6.3) $$

with

$$ S = \lambda_j - \lambda_j^* - (\theta_j - \theta_j^*), \qquad (6.4) $$

where $\theta_j$ is the angle of rotation of body $j$, $\Theta_j' = C_j \omega_j'$ is the conjugated momentum of $\theta_j$, $C_j = \alpha_j m_j R_j^2$ is the moment of inertia of body $j$, $\omega_j' = d\theta_j/dt$ is its rotation rate, and $\alpha_j$ is a dimensionless structure constant depending on the state equation of body $j$ ($\alpha_j = 2/5$ for an homogeneous body). For any quantity $z$ appearing in $U_t^{(j)}$, we denote

$$ z^* = z(t - \Delta t_j), \qquad (6.5) $$

while for any angular quantity $\zeta$ appearing in $U_t^{(j)}$, we define

$$ \Delta\zeta = \zeta - \zeta^* = \zeta(t) - \zeta(t - \Delta t_j). \qquad (6.6) $$

We write the tidal Hamiltonian $H_t$ in the variables $(J, J_2, X_j; \xi, \xi_2, \bar{X}_j)$ introduced in Sect. 3.1.2 and we define $\Theta_j = \Theta_j'/(m\bar{a}^2\eta)$ in order to normalize the action variable $\Theta_j'$ (see Eq. (3.18)). We showed in Sect. 3.1.1 that it is enough to expand the Hamiltonian $\iota H_P$ at order 0 (Eq. (3.11)) in the vicinity of the Keplerian resonance given by Eq. (3.6). However, we show after giving the equations of motions with tides (Eqs. (6.15)), that expanding $U_t^{(j)}$ at order 0 in the neighbourhood of the Keplerian resonance loses almost all tidally relevant dynamics, and we have to keep a dependency on the variables $J$ and $J_2$ in $U_t^{(j)}$. We thus keep exact expressions in these variables. In order to only consider the secular (i.e. long-term) dynamics, we average $U_t^{(j)}$ over $\lambda_j$ and $U_t^{(j)}$ depends on this angle through $\Delta\lambda_j$ only. Similarly, $U_t^{(j)}$ depends on $\theta_j$ through $\Delta\theta_j$ only, in virtue of Eq. (6.4). Since $\iota H_P$ was expanded up to the fourth degree in eccentricity, we also expand
$U_t^{(j)}$ up to that order. The expansion is achieved using the results of Laskar and Robutel (1995), viewing $\mathbf{r}_j$ and $\mathbf{r}_j^*$ as two distinct bodies. We define the quantities

$$q_j = \kappa_2^{(j)} \varphi_j^5, \quad \text{and} \quad \varphi_j = \frac{R_j}{\bar{a}}, \qquad (6.7)$$

and performing the rescaling (3.22), we get for $\mathcal{H}_t = H_t/(m\bar{a}^2\eta)$ the expression<fnref n="1" />

$$\begin{align}
\mathcal{U}_t^{(j)} &= -q_j \eta \frac{m_0}{m} \mathcal{R}_j^{-6} \mathcal{R}_j^{*-6} \left\{ A_t^{(j)} + \mathcal{D}_2^{(j)} + \mathcal{D}_4^{(j)} \right\} = \frac{U_t^{(j)}}{m \bar{a}^2 \eta}, \tag{6.8} \\
\mathcal{T}_j &= \frac{\eta}{2 \alpha_j} \frac{m}{m_j} \frac{\Theta_j^2}{\varphi_j^2} = \frac{T_j}{m \bar{a}^2 \eta}, \nonumber
\end{align}$$

where

$$\mathcal{D}_2^{(j)} = B_t^{(j)} (\mathcal{R}_j^{-1} X_j \bar{X}_j + \mathcal{R}_j^{*-1} X_j^* \bar{X}_j^*) + (\mathcal{R}_j \mathcal{R}_j^*)^{-1/2} (C_t^{(j)} X_j \bar{X}_j^* + \bar{C}_t^{(j)} X_j^* \bar{X}_j), \quad (6.9)$$

and

$$\begin{align}
A_t^{(j)} &= \frac{1}{4} + \frac{3}{4} \cos [2 (\Delta\lambda_j - \Delta\theta_j)], \nonumber \\
B_t^{(j)} &= \frac{3}{8} - \frac{15}{8} \cos [2 (\Delta\lambda_j - \Delta\theta_j)], \tag{6.10} \\
C_t^{(j)} &= \frac{3}{32} e^{i(\Delta\lambda_j - 2\Delta\theta_j)} + \frac{9}{16} e^{-i\Delta\lambda_j} + \frac{147}{32} e^{-i(3\Delta\lambda_j - 2\Delta\theta_j)}. \nonumber
\end{align}$$

The term $\mathcal{D}_4^{(j)}$ contains the fourth order in eccentricity and is given in Appendix B.3. The pseudo-Hamiltonian that we consider for the model is then

$$\mathcal{H} = \mathcal{H}_K(J, J_2) + \iota \mathcal{H}_P(\xi, X_1, X_2, \bar{X}_1, \bar{X}_2) + \mathcal{H}_t, \qquad (6.11)$$

where

$$\mathcal{H}_K = -\frac{3}{2}\eta \frac{m_1 + m_2}{m} J^2 - \frac{3}{2}\eta \frac{m}{m_1 + m_2} J_2^2 + \eta J_2, \qquad (6.12)$$

the perturbation $\iota \mathcal{H}_P = \mathcal{H}^{(0)} + \mathcal{H}^{(2)} + \mathcal{H}^{(4)}$ is given by Eq. (3.29), and

$$\mathcal{H}_t = \sum_{j \in \{1, 2\}} \mathcal{U}_t^{(j)}(J, J_2, J^*, J_2^*, \Delta\xi, \Delta\xi_2, \Delta\theta_j, X_j, X_j^*, \bar{X}_j, \bar{X}_j^*) + \sum_{j \in \{1, 2\}} \mathcal{T}_j(\Theta_j). \quad (6.13)$$

## 6.1.2 Equations of motions in the presence of tides

The equations of motion are derived from the pseudo-Hamiltonian (6.11), following the procedure described in Sect. 5.3 and using Eq. (3.23). Denoting

$$\vartheta_j = 1 - \frac{\omega_j'}{\eta}, \qquad (6.14)$$

<footnote n="1">¹$\mathcal{R}_j(J, J_2)$ is defined by Eq. (3.20) and $m = \sqrt{m_1 m_2}$.</footnote>
The differential system describing the motion of two planar co-orbital planets, undergoing tides, reads<fnref n="2" />

$$
\begin{align}
\dot{\vartheta}_j &= -3\eta\alpha_j^{-1}\frac{m_0}{m_j}\vartheta_j^{-2}\frac{q_j}{Q_j}\mathcal{R}_j^{-12} \{\vartheta_j + 3(1-\mathcal{R}_j) + h_2^{(j)}\mathcal{R}_j^{-1}X_j\bar{X}_j + h_4^{(j)}\mathcal{R}_j^{-2}X_j^2\bar{X}_j^2\}, \nonumber \\
\dot{J} &= -\frac{\partial\mathcal{H}_P}{\partial\xi} + (1-\delta)\dot{J}_2^{(1)} - \delta\dot{J}_2^{(2)}, \nonumber \\
\dot{J}_2 &= \dot{J}_2^{(1)} + \dot{J}_2^{(2)}, \tag{6.15} \\
\dot{\xi} &= \frac{\partial\mathcal{H}_K}{\partial J} + 6\eta q_1 \frac{m_0}{m_1}\mathcal{R}_1^{-13}\mathfrak{Q}_2(\mathcal{R}_1^{-1}X_1\bar{X}_1) - 6\eta q_2 \frac{m_0}{m_2}\mathcal{R}_2^{-13}\mathfrak{Q}_2(\mathcal{R}_2^{-1}X_2\bar{X}_2), \nonumber \\
\dot{X}_j &= -2i\frac{m}{m_j}\frac{\partial(\mathcal{H}^{(2)} + \mathcal{H}^{(4)})}{\partial\bar{X}_j} - 3\eta\frac{q_j}{Q_j}\frac{m_0}{m_j}\mathcal{R}_j^{-13}X_j \left\{ p_2^{(j)} - \frac{5}{2}iQ_j + \frac{X_j\bar{X}_j}{\mathcal{R}_j}(p_4^{(j)} - \frac{65}{4}iQ_j) \right\}, \nonumber
\end{align}
$$

where

$$
\begin{align}
\dot{J}_2^{(j)} &= -3\eta \frac{q_j}{Q_j} \frac{m_0}{m} \mathcal{R}_j^{-12} \{ \vartheta_j + 3(1-\mathcal{R}_j) + k_2^{(j)} \mathcal{R}_j^{-1} X_j \bar{X}_j + k_4^{(j)} \mathcal{R}_j^{-2} X_j^2 \bar{X}_j^2 \}, \nonumber \\
h_2^{(j)} &= \frac{93}{2} + \frac{15}{2}\vartheta_j - \frac{81}{2}\mathcal{R}_j, \qquad h_4^{(j)} = \frac{1989}{8} + \frac{195}{8}\vartheta_j - \frac{819}{4}\mathcal{R}_j, \nonumber \\
k_2^{(j)} &= \frac{157}{2} + \frac{27}{2}\vartheta_j - 69\mathcal{R}_j, \qquad k_4^{(j)} = \frac{2515}{4} + \frac{273}{4}\vartheta_j - \frac{2091}{4}\mathcal{R}_j, \tag{6.16} \\
p_2^{(j)} &= 32 + 6\vartheta_j - \frac{57}{2}\mathcal{R}_j, \qquad p_4^{(j)} = \frac{3041}{8} + \frac{351}{8}\vartheta_j - 318\mathcal{R}_j, \nonumber \\
\mathfrak{Q}_2(Z) &= 1 + \frac{65}{8}Z + \frac{455}{16}Z^2. \nonumber
\end{align}
$$

The set of Eqs. (6.15) does not depend on the angle $\xi_2$. Since we are not interested in the dynamics of this angle, we discarded the line $\xi_2$ from the differential system. If the Hamiltonian $\mathcal{H}_t$ had been expanded at order 0 in the vicinity of the Keplerian resonance (3.6) (i.e. at order 0 in the variables $J$ and $J_2$), then the subsequent differential system would be the set of Eqs. (6.15), with the substitution $\mathcal{R}_j = 1$. Let us study the circular dynamics in that case, that we treated in the point-mass approximation in Sect. 3.2.1. We show in Sect. 6.2.2 that the rotation rates $\omega_j$ are damped by tides in a timescale much smaller than the other tidal timescales, so we also perform the substitution $\vartheta_j = 0$ in Eqs. (6.15) to consider the dynamics after reaching the spin-orbit synchronization. In the circular case, we substitute $X_j$ by 0 and we obtain, for the librating angle $\xi$, the differential equation

$$
\ddot{\xi} = -3\eta^2 \iota \sin \xi (1 - (2 - 2 \cos \xi)^{-3/2}). \quad (6.17)
$$

All tidal contributions have disappeared and we end up with the differential equation describing the motion of $\xi$ in the point-mass approximation (Yoder et al., 1983, Eq. (9)). This equation predicts a periodic motion for $\xi$ and no long-term effects due to tides are described. It is thus necessary to keep a dependency in $J$ and $J_2$ in the tidal Hamiltonian to have a satisfying description of the tidal dynamics of the co-orbital motion.

<footnote n="2">²$\delta$ was defined by Eq. (3.14).</footnote>
### 6.1.3 Conservation of the total angular momentum with tides

While tides do not preserve the total energy of the system because of the dissipation, the total angular momentum of the system star–planet–planet has to be conserved, even with tides, since the system is closed. We can show that the differential system (6.15) is indeed consistent with the conservation of the total angular momentum. The normalized<fnref n="3" /> total angular momentum $\mathcal{C}$ reads (Robutel and Pousse, 2013, Eq. (6))

$$ \mathcal{C} = \sum_j \frac{m_j}{m} \mathcal{R}_j - \frac{1}{2} \sum_j \frac{m_j}{m} X_j \bar{X}_j + \sum_j \alpha_j \frac{m_j}{m} \varphi_j^2 (1 - \vartheta_j). \quad (6.18) $$

From (6.15) we get

$$ \dot{\mathcal{C}} = \sum_{j=1}^{2} 3\eta \frac{q_j}{Q_j} \frac{m_0}{m} \mathcal{R}_j^{-13} X_j \bar{X}_j \left\{ h_2^{(j)} - k_2^{(j)} + p_2^{(j)} + \mathcal{R}_j^{-1} X_j \bar{X}_j (h_4^{(j)} - k_4^{(j)} + p_4^{(j)}) \right\}, \quad (6.19) $$

and the total angular momentum is conserved since we have (Eq. (6.16))

$$ \begin{gathered} h_2^{(j)} - k_2^{(j)} + p_2^{(j)} = 0 \quad \text{and} \\ h_4^{(j)} - k_4^{(j)} + p_4^{(j)} = 0. \end{gathered} \quad (6.20) $$

## 6.2 Equilibria and linearization

In this section, we compute the equilibria of the differential system (6.15) and the eigenvalues in their vicinity, in order to infer the stability of the system, using the results of Sect. 2.4.1.

### 6.2.1 The eigenvalues

Since the system (6.15) is a perturbation of the Hamiltonian system derived from $\mathcal{H}_K+\iota\mathcal{H}_P$, its fixed points are a perturbation of the Lagrangian equilibria. Only the tidal contributions of first order in $\iota^{-1}(q_1+q_2)$ appear in (6.15), and it is enough to restrict to the first order to compute the position of the fixed points. We find, for $j \in \{1, 2\}$,

$$ \begin{gathered} \vartheta_j = 0, \\ \xi - \frac{\pi}{3} = 0, \\ -3 \frac{m_1 + m_2}{m} J + 6q_1 \frac{m_0}{m_1} - 6q_2 \frac{m_0}{m_2} = 0, \\ X_j = 0. \end{gathered} \quad (6.21) $$

Note that this choice for the fixed points does not make the right hand side of system (6.15) to be exactly zero, but rather a quantity of second order in $\iota^{-1}(q_1+q_2)$. Infinitely many choices are possible for Eq. (6.21) that make the right-hand side of (6.15) a quantity of second order in $\iota^{-1}(q_1+q_2)$, and we choose one of them. This choice is the most natural

<footnote n="3">³Normalized by $m\bar{a}^2\eta$.</footnote>
because it guarantees that the fixed points correspond to a solid rotation of the whole system ($\vartheta_j = 0$). Equations (6.21) only provide 6 independent equations for 7 variables<fnref n="4" />. In order to have a unique fixed point in the neighbourhood of which to linearize the system, we arbitrarily choose

$$f_1 + f_2 = 0, \qquad (6.22)$$

where the $f_j$ are defined by Eq. (3.21). At first order in $\iota^{-1}(q_1 + q_2)$, the linearized system does not depend on the choice of Eqs. (6.21) & (6.22).

Before we linearize the differential system, we show that our model allows us to recover a well known result, called pseudo-synchronization. We fix to 0 the value of $\dot{\vartheta}_j$ in (6.15). If we place ourselves at the Keplerian resonance; that is, at $a_j = \bar{a}$ (or equivalently $\mathcal{R}_j = 1$), we find<fnref n="5" /> for the equilibrium value of the rotation rate of planet j

$$\frac{\omega'_j}{\eta} = 1 + 6e_j^2 + \frac{3}{8}e_j^4 + \mathcal{O}(e_j^6). \qquad (6.23)$$

This is the pseudo-synchronization of the rotation rates. At non-zero eccentricities, the rotation rates $\omega'_j$ are not synchronized with $\dot{\lambda}_j$, but rather slightly super-synchronous. The same result can be reached by equating to zero the right-hand side of Eq. (3) in Correia and Laskar (2004). Equation (6.23) and Sect. 6.1.3 provide two independent verifications of the set of Eqs. (6.15).

Let $\mathcal{X}_0 = {}^t(\vartheta_{1,0}, \vartheta_{2,0}, \xi_0, J_0, J_{2,0}, X_{1,0}, X_{2,0})$ be the unique solution of Eqs. (6.21) & (6.22), that is, the equilibrium value of $\mathcal{X} = {}^t(\vartheta_1, \vartheta_2, \xi, J, J_2, X_1, X_2)$. Once linearized in the vicinity of $\mathcal{X}_0$, the differential system (6.15) reads

$$\frac{d}{dt} \Delta \mathcal{X} = (Q_0 + Q_1) \Delta \mathcal{X}, \qquad (6.24)$$

where $\Delta \mathcal{X} = \mathcal{X} - \mathcal{X}_0$. The matrix $Q_0$ derives from the conservative Hamiltonian $\mathcal{H}_K + \iota \mathcal{H}_P$, while $Q_1$ corresponds to the tidal contribution $\mathcal{H}_t$. The $7 \times 7$ matrix $Q_0 + Q_1$ is block diagonal with a $5 \times 5$ block corresponding to $(\vartheta_1, \vartheta_2, \xi, J, J_2)$, and $2 \times 2$ block, corresponding to $X_1$ and $X_2$. This shows that, even in the presence of tides, the circular and eccentric dynamics are uncoupled near the Lagrangian equilibria $L_{4,5}$. A detailed expression of these matrices is given in appendix C.3. The set of eigenvalues of $Q_0$ is

$$\{0, 0, 0, i\nu, -i\nu, ig_1, ig_2\}, \qquad (6.25)$$

where (Eqs. (3.32) & (3.42))

$$\nu = \eta \sqrt{\frac{27\iota}{4}}, \quad g_1 = \eta \frac{27\iota}{8}, \quad g_2 = 0. \qquad (6.26)$$

Among the first three 0 eigenvalues, two correspond to the constant rotation rate of the planets and another to the conservation of the total angular momentum (constancy of $J_2$). The $\pm i\nu$ eigenvalues give the libration frequency of the resonant angle $\xi = \lambda_1 - \lambda_2$ around the Lagrangian equilibria $L_{4,5}$ (Sect. 3.2.1), and the last two eigenvalues give the precession frequency of the pericentres in the eigenmodes anti-Lagrange for $g_1$ and Lagrange for $g_2$

<footnote n="4">⁴This is due to the conservation of the total angular momentum of the system.</footnote>

<footnote n="5">⁵Using $X_j \bar{X}_j = 2(1 - \sqrt{1 - e_j^2}) = e_j^2 + e_j^4/4 + \mathcal{O}(e_j^6)$.</footnote>
(Sect. 3.2.2). In particular, all seven eigenvalues of $Q_0$ are pure imaginary and without the contribution of tides, the trajectories of the linearized system are quasiperiodic.

The matrix $Q_0 + Q_1$ is too gruesome for its eigenvalues to be computed by finding the roots of its characteristic polynomial. Instead, since $Q_1$ is a small perturbation of $Q_0$, and because we can easily find a diagonal basis of $Q_0$, we use the perturbative approach described in Sect. 2.4.2 to find the eigenvalues of $Q_0 + Q_1$ at first order in $\iota^{-1} (q_1 + q_2)$. The computations are performed with the algebraic manipulator Maxima and the set of eigenvalues of $Q_0 + Q_1$ is

$$ \{\lambda_1, \lambda_2, 0, \lambda, \bar{\lambda}, \lambda_{\text{AL}}, \lambda_{\text{L}}\}, \qquad (6.27) $$

with

$$
\begin{align}
\lambda_j &= \left[ -3\eta\alpha_j^{-1}\frac{q_j}{Q_j}\varrho_j^{-2}\frac{m_0}{m_j} + 9\eta\iota^{-1}\frac{q_j}{Q_j} \right] < 0, \nonumber \\
\lambda &= \left[ \frac{9}{2}\eta\iota^{-1}\left(\frac{m_1}{m_2}\frac{q_2}{Q_2} + \frac{m_2}{m_1}\frac{q_1}{Q_1}\right) + i\nu \left[ 1 + 13\iota^{-1}\left(\frac{m_1}{m_2}q_2 + \frac{m_2}{m_1}q_1\right) \right] \right], \nonumber \\
\lambda_{\text{AL}} &= \left[ -\frac{21}{2}\eta\iota^{-1}\left(\frac{m_1}{m_2}\frac{q_2}{Q_2} + \frac{m_2}{m_1}\frac{q_1}{Q_1}\right) + ig_1 \left[ 1 + \frac{20}{9}\iota^{-2}\left(\frac{m_1}{m_2}q_2 + \frac{m_2}{m_1}q_1\right) \right] \right], \tag{6.28} \\
\lambda_{\text{L}} &= \left[ -\frac{21}{2}\eta\iota^{-1}\left(\frac{q_1}{Q_1} + \frac{q_2}{Q_2}\right) + \frac{15}{2}i\eta\iota^{-1}(q_1 + q_2) \right]. \nonumber
\end{align}
$$

We note that the eigenvalues are no longer pure imaginary and we boxed the real parts for a better visualization. The real parts are proportional to the inverse quality factors $Q_j^{-1}$, while the perturbations to the imaginary parts are not. As a consequence, elastic tides, for which $Q_j = \infty$ (or equivalently, $\Delta t_j = 0$, see Eq. (5.44)), do not yield dissipation but only change slightly the fundamental frequencies of the system. We computed the eigenvectors of the linearized system and we show in appendix C.3 how the eigenmodes Lagrange and anti-Lagrange are slightly modified by tides.

The matrix $Q_0 + Q_1$ is singular, and one of its eigenvalues is zero, corresponding to the conservation of the total angular momentum (Sect. 6.1.3). As a result of tidal dissipation, the eigenvalues $\lambda_{\text{AL}}$ and $\lambda_{\text{L}}$, perturbations of $ig_1$ and $ig_2$, respectively, have non-zero negative real parts. Therefore, both eccentric eigenmodes Lagrange and anti-Lagrange are damped to zero. Similarly, $\lambda_j < 0$, and the rotation rates of both planets are also damped. On the contrary, the real part of $\lambda$ and $\bar{\lambda}$, perturbations of $i\nu$ and $-i\nu$, are strictly positive, which leads to an exponential increase of the libration amplitude of $\xi$ when the system is close to the Lagrangian equilibria $L_4$ or $L_5$.

The fact that $\Re \lambda$ is strictly positive is not a surprise. Indeed, $L_4$ and $L_5$ are maxima of energy in the phase space (Fig. 3.1), and since tides are responsible for a decrease in the orbital energy, we expect the system to be driven away from the Lagrangian equilibria by tides, hence the positivity of $\Re \lambda$. The fact that $L_4$ and $L_5$ are maxima of energy is not a surprise either, as Moeckel (2017) proved that a system of $N \ge 3$ rigid bodies rotating rigidly and uniformly around their axis of total angular momentum cannot be an energy minimizer of the space phase. Therefore, we could have expected without calculation that the elliptic Lagrangian equilibria are energy maximizers$^6$, and as a consequence, that

<footnote>$^6$They must be critical points, in order to be equilibria of the averaged model, and they cannot be</footnote>
dissipative tides make them unstable.

## 6.2.2 The timescales of tidal evolution

We define here the characteristic timescale of a given proper mode of the system (6.24) as the time needed for its amplitude to be multiplied (or divided, if the corresponding real part is negative) by a factor $\exp(1)$. According to the eigenvalues (6.28), these times are

$$ 
\begin{aligned}
\tau_{\text{rot}}^{(j)} &= \frac{1}{6\pi} \alpha_j \varphi_j^2 \frac{m_j Q_j}{m_0 q_j} \left[ 1 + \frac{3m_j}{m_1 + m_2} \varphi_j^2 \alpha_j \right] T, \\
\tau_{\text{L}} &= \frac{\iota}{21\pi} \left( \frac{q_1}{Q_1} + \frac{q_2}{Q_2} \right)^{-1} T, \\
\tau_{\text{AL}} &= \frac{\iota}{21\pi} \left( \frac{m_2 q_1}{m_1 Q_1} + \frac{m_1 q_2}{m_2 Q_2} \right)^{-1} T, \\
\tau_{\text{lib}} &= \frac{7}{3} \tau_{\text{AL}},
\end{aligned}
\quad (6.29) $$

where $T = 2\pi/\eta$ is the orbital period. We note that the times $\tau_{\text{rot}}^{(j)}$ are much smaller than the three other characteristic timescales, due to the presence of the factor $\varphi_j^2 \ll 1$. That is, regardless of the parameters and initial conditions, the rotations of the planets are damped to their equilibrium value (6.23) in a timescale such that the eccentricities and the libration angle do not undergo significant damping or excitation. We know from the analysis of the eigenvalues (6.28) that the two eccentric eigenmodes Lagrange and anti-Lagrange are damped to zero, while the libration amplitude of the resonant angle $\xi$ exponentially increases. We now compare the timescales $\tau_{\text{AL}}$, $\tau_{\text{L}}$ and $\tau_{\text{lib}}$ to determine which proper mode amplitude evolves faster. Even though both eccentric eigenmodes are damped, the damping times may be different. By comparing them, we can find if the system favours the Lagrange or the anti-Lagrange configuration. Indeed if

* $\frac{\tau_{\text{AL}}}{\tau_{\text{L}}} < 1$, then the system settles in Lagrange, whereas if

* $\frac{\tau_{\text{AL}}}{\tau_{\text{L}}} > 1$, then the system settles in anti-Lagrange.

Moreover, comparing the time $\tau_{\text{lib}}$ with the eccentric times $\tau_{\text{AL}}$ and $\tau_{\text{L}}$ allows us to determine if the system is still eccentric or already nearly circular when the libration amplitude has significantly increased. We have

$$ \tau_{\text{lib}} = \frac{\iota}{9\pi\Omega} \frac{\tau_{\text{AL}}}{\tau_{\text{L}}} T, \quad (6.30) $$

where we defined $\Omega = q_1/Q_1 + q_2/Q_2$ the sum of the dissipation rates. Equation (6.30) shows that, for a given sum of the planetary masses, $\iota$, and sum of the dissipation rates, $\Omega$, the system moves away from $L_{4,5}$ faster if it favours Lagrange, and slower if it favours anti-Lagrange.

These results are a priori valid only in a small neighboorhoud of $L_{4,5}$, but in fact, the simulations from Sect. 6.3 show that the behaviour of the system near $L_{4,5}$ is valid

<footnote>saddle points, or they would be hyperbolic.</footnote>
even at high libration amplitudes. This means that the exponential increase in libration amplitude is unbounded. Once the libration amplitude of $\xi$ increased too much, Eq. (3.37) is not verified anymore, the co-orbitals enter each other's Hill sphere and planet-planet interactions disrupt the system. The planets either get ejected from the co-orbital resonance, or they collide. The disruption time depends on $\tau_{\text{lib}}$, and so, Lagrange-like systems have a short life expectancy, while anti-Lagrange-like systems have a long life expectancy (for a given $\iota$ and $\Omega$). Moreover, a Lagrange-like system is eccentric when old$^7$ (as long as it was eccentric when young), while an anti-Lagrange-like system is always circular when old. Indeed, for an anti-Lagrange-like system, the characteristic time $\tau_e$ of eccentricity damping is given by $\tau_{\text{AL}}$ and then, $\tau_{\text{lib}} = 7\tau_e/3 \sim \tau_e$ ensures that the eccentricity is damped when the libration amplitude has significantly increased. On the other hand, for a Lagrange-like system, $\tau_e = \tau_L$ and then, $\tau_{\text{lib}} \ll \tau_e$ ensures that the system is still eccentric when the libration amplitude has significantly increased.

Although it is clear that $\tau_{\text{lib}}$ depends on the semi-major axis $\bar{a}$, Eq. (6.30) does not explicitly show it, since the dependency on $\bar{a}$ is hidden in $q_j$, $Q_j$ and $T$. We have

$$ \tau_{\text{lib}} \propto \frac{\tau_{\text{AL}}}{\tau_{\text{L}}} \bar{a}^{\beta}, \qquad (6.31) $$

where $\beta$ is 8 for the constant-$\Delta t$ model and 6.5 for the constant-$Q$ model. Thus, for a large $\bar{a}$ or a strong anti-Lagrange tendency, co-orbital planetary systems may survive for the entire life-time of the star in the main sequence.

Interestingly, the ratio between the two eccentric damping timescales depends only on the ratio between the planetary masses and the ratio between the dissipation rates inside the planets. That is, denoting

$$ x = \frac{m_1}{m_2} \quad \text{and} \quad y = \frac{q_2/Q_2}{q_1/Q_1} = \frac{\kappa_2^{(2)} \varphi_2^5/Q_2}{\kappa_2^{(1)} \varphi_1^5/Q_1}, \qquad (6.32) $$

we have

$$ \frac{\tau_{\text{AL}}}{\tau_{\text{L}}} = \frac{x(1+y)}{1+yx^2}. \qquad (6.33) $$

The equality between the eccentric damping timescales ($\tau_{\text{AL}} = \tau_{\text{L}}$) occurs at $x = 1$ and $xy = 1$, plotted by black lines in Fig. 6.1, where we also show the ratio $\tau_{\text{AL}}/\tau_{\text{L}}$ as a function of $x$ and $y$. We clearly observe two regions, corresponding to Lagrange-like and anti-Lagrange-like systems.

Let us assume that the quality factor $Q_j$ is mass independent. We then deduce from $\varphi_j \propto m_j^{1/3}$ and Eq. (6.32) that

$$ y \propto x^{-5/3}. \qquad (6.34) $$

The blue line in Fig. 6.1 plots $y = x^{-5/3}$. It shows the path followed by a system with two initially identical planets (white spot) when we change the mass repartition between them. We conclude that for a given sum of the planetary masses $\iota$, and sum of the dissipation rates $\Omega$, the expectancy of life of the system is at its shortest for $m_1 = m_2 = \iota m_0/2$ (orange area) and it tends towards infinity if either $m_1$ or $m_2$ tends towards $\iota m_0$, or towards zero (yellow/white area).

<footnote>$^7$Old means that its age is significant with respect to its expectancy of life.</footnote>
<figure><img src="image_14.png" /><figcaption>Fig. 6.1 — Value of $\tau_{AL}/\tau_L$ in colorscale. The black straight lines are the locations of the points where $\tau_{AL} = \tau_L$. The blue line is the path followed by a system of variable mass repartition. The black spots are the positions of the numerical simulations of Sect. 6.3. The white dot ($x = 1, y = 1$) is where Rodríguez et al. (2013) performed all of their simulations. It complies with $\tau_{AL} = \tau_L$, explaining why they did not see any hierarchy between the eccentric eigenmodes. Systems in the red-black regions settle into Lagrange, have a short life expectancy and are eccentric when old, while systems in the yellow-white regions settle into anti-Lagrange, have a long life expectancy and are circular when old.</figcaption></figure>

Combining Eqs. (6.30) & (6.33), the characteristic timescale for the increase of the libration amplitude reads

$$ \tau_{\text{lib}} = \frac{\iota}{9\pi\Omega} \frac{x(1+y)}{1+yx^2} T. \qquad (6.35) $$

In the case of identical co-orbital planets (same mass, radius and tidal parameters), we derive for the characteristic timescale of destruction the simple expression

$$ \tau_{\text{lib}} = \frac{m_j Q_j T}{9\pi m_0 q_j} = 0.565 \text{ Gyr} \left(\frac{m_j}{M_{\oplus}}\right) \left(\frac{M_{\odot}}{m_0}\right)^{3/2} \left(\frac{\bar{a}}{0.04 \text{ AU}}\right)^{6.5} \left(\frac{R_{\oplus}}{R_j}\right)^5 \left(\frac{0.0011}{\kappa_2^{(j)}/Q_j}\right), \quad (6.36) $$

where we considered two co-orbital Earth-like planets for the numerical evaluation in Eq. (6.36). The value of $\kappa_2/Q$ for the Earth is given in Table 5.1.

### 6.2.3 Looking for co-orbital exoplanets

We can use the results from Sect. 6.2.2 to estimate whether or not an already discovered exoplanet may have an undetected co-orbital companion. For a co-orbital system, we denote $\tau_{hs}$ and $\tau_{dest}$ the time needed to reach the horseshoe-shaped orbits (that is, to cross the separatrix emanating from $L_3$), and the time needed for close encounters to disrupt
<figure><img src="image_15.png" /><figcaption>Fig. 6.2 — Disruption time (Gyr) of an hypothetical co-orbital system as a function of the orbital period of the observed planet and of the mass of the host star. The couple Saturn+Earth is plotted in orange while the couple Earth+Earth is plotted in blue. The solid lines plot the minimum of the main sequence duration and of the age of the universe for both couples. Systems below these lines may have been already destroyed at the time of observation, but systems above outlive either their host star or the age of the universe. The blue and orange dots correspond to HD 158259 c and HD 102956 b, respectively. Note that, for a given orbital period, co-orbital systems live longer around a massive host star because they have a larger semi-major axis.</figcaption></figure>

the co-orbital resonance, respectively. In the rest of this section, we use the accessible observational parameters to approximate $\tau_{\text{dest}}$ and then discard cases such that $\tau_{\text{dest}}$ is too small.

The libration amplitude $\Delta\xi = \max(\xi - 60^{\circ})$ is defined as the angular distance to $L_4$. If $\Delta\xi_0$ is the initial libration amplitude of the system and if $\Delta\xi$ is a small libration amplitude greater than $\Delta\xi_0$, then, by definition of $\tau_{\text{lib}}$, the time $\tau_{\xi}$ needed to reach the libration amplitude $\Delta\xi$ is given by

$$ \tau_{\xi} = \ln \left( \frac{\Delta\xi}{\Delta\xi_0} \right) \tau_{\text{lib}}. \qquad (6.37) $$

Equation (6.37) cannot be used, a priori, to predict the time $\tau_{\text{hs}}$, since the separatrix emanating from $L_3$ is far from the fixed point $L_4$. Nevertheless, by performing numerical simulations of the set of Eqs. (6.15) with arbitrary parameters and initial conditions, we verify that the expression

$$ \tau_{\text{hs}} \approx \ln \left( \frac{60^{\circ}}{\Delta\xi_0} \right) \tau_{\text{lib}}, \qquad (6.38) $$

is always a good approximation as long as $\Delta\xi_0 \le 15^{\circ}$, where $\Delta\xi_0$ is in arc degrees. Although $\tau_{\text{hs}}/\tau_{\text{lib}}$ only depends on $\Delta\xi_0$, we expect that $\tau_{\text{dest}}/\tau_{\text{lib}}$ also depends on $\iota$, since this parameter controls the maximum libration amplitude before the system becomes
<table><caption>Table 6.1 — Time to horseshoe-shaped orbits $\tau_{hs}$ for a variety of co-orbital systems, computed with Eq. (6.40). The tidal parameters are those of Table 5.1 and only the six bodies for which $\kappa_2/Q$ is well constrained are included. In general, $\tau_{dest}$ is close to $\tau_{hs}$.</caption><thead><tr><th>Co-orbital pair</th><th>τ<sub>hs</sub> (Gyr)</th><th>Co-orbital pair</th><th>τ<sub>hs</sub> (Gyr)</th><th>Co-orbital pair</th><th>τ<sub>hs</sub> (Gyr)</th></tr></thead><tbody><tr><td>Earth &amp; Earth</td><td>3.612</td><td>Moon &amp; Mars</td><td>28.48</td><td>Mars &amp; Saturn</td><td>5.000</td></tr><tr><td>Earth &amp; Moon</td><td>43.81</td><td>Moon &amp; Jupiter</td><td>50.63</td><td>Jupiter &amp; Jupiter</td><td>0.7214</td></tr><tr><td>Earth &amp; Mars</td><td>5.552</td><td>Moon &amp; Io</td><td>4.374</td><td>Jupiter &amp; Io</td><td>2.072</td></tr><tr><td>Earth &amp; Jupiter</td><td>3.567</td><td>Moon &amp; Saturn</td><td>43.13</td><td>Jupiter &amp; Saturn</td><td>0.04741</td></tr><tr><td>Earth &amp; Io</td><td>2.085</td><td>Mars &amp; Mars</td><td>5.892</td><td>Io &amp; Io</td><td>2.072</td></tr><tr><td>Earth &amp; Saturn</td><td>1.803</td><td>Mars &amp; Jupiter</td><td>5.878</td><td>Io &amp; Saturn</td><td>2.054</td></tr><tr><td>Moon &amp; Moon</td><td>50.76</td><td>Mars &amp; Io</td><td>2.250</td><td>Saturn &amp; Saturn</td><td>0.03704</td></tr></tbody></table>

unstable (Eq. (3.37)). The smaller $\iota$ is, the larger the maximum libration amplitude is. Since $\iota = 3 \times 10^{-4}$ is the highest value allowing horseshoe-shaped orbits (Leleu et al., 2015, Fig. 4), the equality $\tau_{dest} = \tau_{hs}$ occurs at this value. If we take $\Delta\xi_0 = 10^\circ$, we have

$$1/2 \tau_{hs} \le \tau_{dest} \le 2 \tau_{hs} \Leftrightarrow 10^{-9} \lesssim \iota \lesssim 0.005, \qquad (6.39)$$

and so, $\tau_{hs}$ and $\tau_{dest}$ do not differ by more than a factor 2 for a wide range of $\iota$. We thus consider $\tau_{dest} \approx \tau_{hs}$ in this range and Eqs. (6.35) & (6.38) can be used to predict $\tau_{dest}$.

Figure 6.2 shows the disruption times of a Saturn-Earth co-orbital system (in orange), and of a system of two Earth-like planets (in blue), with $\Delta\xi_0 = 10^\circ$. The tidal parameters adopted for these systems are given in Table 5.1. The solid line plots min($\tau_{ms}$, $\tau_u$) where $\tau_{ms} = 10^{10}$ yr ($m_0/M_\odot)^{-2}$ is the duration of the main-sequence of the host star and $\tau_u = 13.77$ Gyr is the age of the Universe. Assuming that the Earth and Saturn are representative of the average rocky planet and gas giant, Fig. 6.2 tells us if an already detected gas giant may have a companion (orange lines), or if an already detected rocky planet may have a companion (blue lines). The detected exoplanet is located in Fig. 6.2 using its orbital period and the mass of the host star. If it is below its associated solid line (orange for a gas giant, blue for a rocky planet), then its companion was already ejected, if it ever existed. On the contrary, if it is above this line, the hypothetical pair of co-orbital outlives either the host star or the age of the Universe and it is worth looking for the companion. As an example, it is very unlikely to find a co-orbital companion for the rocky planet HD 158259 c (Hara et al., 2020), plotted with a blue dot in Fig. 6.2, but we cannot rule out that the gas giant HD 102956 b (Luhn et al., 2019), plotted with an orange dot in Fig. 6.2, has a co-orbital companion.

Combining Eqs. (6.35) & (6.38), the time $\tau_{hs}$ to reach the horseshoe-shaped orbits takes the form

$$\tau_{hs} = \frac{\iota x (1 + y)}{9\pi\Omega (1 + yx^2)} \ln \left( \frac{60^\circ}{\Delta\xi_0} \right) T. \qquad (6.40)$$

We can use the tidal parameters given by Table 5.1 to compute this time for a large variety of co-orbital systems. We give in Table 6.1 the time $\tau_{hs}$ for hypothetical co-orbital pairs of exoplanets made up of Solar System bodies. The semimajor axis is $\bar{a} = 0.04$ AU and the mass of the host star is $m_0 = M_\odot$, but $\tau_{hs}$ can easily be deduced for other values using the exponents of Eq. (6.36), that are still valid for different co-orbitals. We choose
$\Delta\xi_0 = 0.1^\circ$, and again, it is straightforward to extend the results to another choice of $\Delta\xi_0$ using Eq. (6.38).

## 6.3 Numerical simulations

In this section, in order to verify the analytical results of Sect. 6.2, we perform some numerical simulations of planetary systems representative of the different dynamical regimes of co-orbital planets undergoing tidal interactions with the star, such as those described in Sect. 6.2.

### 6.3.1 Procedure

We numerically integrate the six systems that correspond to the black dots on Fig. 6.1. For all of them we choose the initial conditions and parameters<fnref n="8" />

$$ 
\begin{array}{l@{\quad}l@{\quad}l}
\iota = 2 \times 10^{-4}, & \Omega = 4 \times 10^{-13}, & m_0 = M_{\odot}, \\
\rho_2 = 4\rho_1 = 2000 \, \text{kg/m}^3, & e_{2,0} = 2e_{1,0} = 0.04, & \alpha_{1,0} = \alpha_{2,0} = 0.33, \\
a_{1,0} = a_{2,0} = \bar{a} = 0.02 \, \text{AU}, & \vartheta_{1,0} = \vartheta_{2,0} = 0, & \varpi_{1,0} = \varpi_{2,0} = 0, \quad (6.41) \\
& & \xi_0 = 62^\circ,
\end{array}
$$

that is, the systems are initially 2° away from $L_4$. Since the total planetary masses $\iota$ and the total dissipation rate $\Omega = q_1/Q_1 + q_2/Q_2$ are the same for all systems, their positions can all be plotted in Fig. 6.1 and their tidal timescales are entirely determined by the values of the mass ratio $x$ and the dissipation rate ratio $y$ (see Eq. (6.32)), which are the only variable parameters between the systems (see Table 6.2). According to Eq. (6.30), we expect that the systems that are in the yellow regions of Fig. 6.1 live longer than those in the red regions. This choice for the initial eccentricities and longitude of the pericentres guarantees that the systems are not initially collinear to the Lagrange or anti-Lagrange configuration, and their evolution allows us to determine if they are Lagrange-like or anti-Lagrange-like. The timescales of the six systems, in number of orbital periods and deduced from Sect. 6.2.2, are given in Table 6.2. The time $\tau_{\text{lib}}$ is the characteristic timescale of libration amplitude excitation given by Eqs. (6.35), while $\tau_{\text{lib}}^{(\text{num})}$ is directly computed from the numerical value of $Z_0 + Z_1$ (given in appendix C.3), without using the perturbative method of Sect. 2.4.2.

<table><caption>Table 6.2 — Values of $x$, $y$ and of the corresponding timescales of the systems. The timescales are given in number of orbital periods, that is, in units of $2\pi/\eta$.</caption><thead><tr><th>#</th><th>color</th><th>x</th><th>y</th><th>τ<sub>lib</sub></th><th>τ<sub>lib</sub><sup>(num)</sup></th><th>τ<sub>AL</sub></th><th>τ<sub>L</sub></th></tr></thead><tbody><tr><td>1</td><td>blue</td><td>10</td><td>100</td><td>1 785 893</td><td>1 845 021</td><td>765 382</td><td>7 578 807</td></tr><tr><td>2</td><td>green</td><td>1/500</td><td>100</td><td>3 570 716</td><td>4 198 710</td><td>1 530 306</td><td>7 578 807</td></tr><tr><td>3</td><td>red</td><td>100</td><td>1/50</td><td>8 973 910</td><td>9 161 859</td><td>3 845 961</td><td>7 578 807</td></tr><tr><td>4</td><td>purple</td><td>100</td><td>1/200</td><td>34 847 651</td><td>34 893 952</td><td>14 934 707</td><td>7 578 807</td></tr><tr><td>5</td><td>yellow</td><td>1/10</td><td>100</td><td>89 303 607</td><td>89 304 263</td><td>38 272 974</td><td>7 578 807</td></tr><tr><td>6</td><td>black</td><td>100</td><td>10<sup>-5</sup></td><td>1 607 641 764</td><td>1 607 642 323</td><td>688 989 327</td><td>7 578 807</td></tr></tbody></table>

<footnote n="8">⁸The variable $\rho_j$ is the density of planet $j$.</footnote>
<figure><img src="image_16.png" /><figcaption>Fig. 6.3 — Value of the libration angle $\xi = \lambda_1 - \lambda_2$ against time, integrated with Eqs. (6.15) (top) and Eqs. (C.13) (bottom). The switching from tadpole to horseshoe-shaped orbits (crossing of the separatrix emanating from $L_3$) is very visible, as the libration amplitude suddenly increases, giving this boot-shape to the plots. As expected from Fig. 6.1 and Table 6.2, systems 4, 5 and 6 live longer than systems 1, 2 and 3.</figcaption></figure>

For some systems, especially system 2, there is a slight difference between $\tau_{\text{lib}}^{(\text{num})}$ and $\tau_{\text{lib}}$. For these systems, the values of $a_1$ and $a_2$ in the matrix $Z_1$ (appendix C.3) are much larger than the other entries of the matrix, and are such that $a_j \sim \nu$, which means that the condition $|\lambda_i - \lambda_j| \gg \epsilon$ in Eq. (2.88) of Sect. 2.4.2 is poorly respected. A smaller value for the sum of the dissipation rates $\Omega$ provides a perfect agreement between $\tau_{\text{lib}}^{(\text{num})}$ and $\tau_{\text{lib}}$ for all six systems, but it also leads to much longer simulations. Here, we purposefully choose a large value for $\Omega$ in order to have quick simulations. On the other hand, there is no disagreement between the analytical values of $\tau_{\text{AL}}$ and $\tau_{\text{L}}$ given by Eq. (6.29) and their numerical counterparts, since the eccentricities are uncoupled from the rest of the variables in the linearized system and the entries of $M_1$ in appendix C.3 respect the condition of Eq. (2.88) for all six systems.

Each system is numerically integrated using two different sets of equations. In the first set, we use the secular Eqs. (6.15) derived in Sect. 6.1.2. In the second set, we use a $n$-body unaveraged and unexpanded model in the regular position-speed coordinates, with the constant-$\Delta t$ model. The second set is given by Eq. (C.13) of appendix C.4. The results are displayed in Figs. 6.3, 6.4 and 6.5. In each figure, the top panels correspond to the set of Eqs. (6.15) while the bottom panels correspond to the set of Eqs. (C.13). All systems except system 6 are integrated long enough for the co-orbital configuration to be destroyed.
<figure><img src="image_17.png" /><figcaption>Fig. 6.4 — Value of $\varpi_1 - \varpi_2$ against time, integrated with Eqs. (6.15) (top) and Eqs. (C.13) (bottom). Blue, green and red systems settle into the Lagrange configuration and have a short life while purple, yellow and black systems settle into the anti-Lagrange configuration and live long. This is in agreement with what we deduce from the analytical results, especially Fig. 6.1 and Table 6.2.</figcaption></figure>

### 6.3.2 Comparison between the model and the complete system

We observe that there is always a very good agreement between Eqs. (6.15) and Eqs. (C.13), except when the libration amplitude is near 360°. Indeed, as the libration amplitude increases, close encounters between the planets mean that planet–planet interactions are no longer perturbations of the Keplerian motion and Eq. (3.37) is not verified anymore. On the top plot of Fig. 6.3, the amplitude of $\xi$ tends towards 360° as time goes to infinity but never reaches it, while on the bottom plot, there exists a finite time when $\xi$ reaches 360°, meaning it is a circulating angle and the planets either collide or get ejected off the co-orbital resonance. All simulations confirm the destruction of the co-orbital resonance. We also note that for small eccentricities, the averaged Eqs. (6.15) (plotted on top) differ from the complete Eqs. (C.13) (plotted at the bottom). We believe that this is due to short period influences that were averaged out in our model. Indeed, while the manifold ($X_1 = X_2 = 0$) is stable by the flow of the averaged Hamiltonian, it is not stable by the flow of the complete Hamiltonian. This may explain the differences between both plots at low eccentricities. This does not discredit though the theoretical results obtained in Sect. 6.2.

Let $\tau_{hs}^{(a)}$, $\tau_{hs}^{(s)}$ and $\tau_{hs}^{(c)}$ be the times needed to reach the horseshoe-shaped orbits according to the analytical expression<fnref n="9" /> (6.38), the simulation of the secular Eqs. (6.15) and the

<footnote n="9">⁹Due to tides being too strong ($\Omega$ is too large) for systems 1 and 2, we considered $\tau_{\text{lib}}^{(\text{num})}$ instead of $\tau_{\text{lib}}$. With weaker tides, there is no difference between these two timescales.</footnote>
<figure><img src="image_18.png" /><figcaption>Fig. 6.5 — Eccentricities of both planets against time, integrated with Eqs. (6.15) (top) and Eqs. (C.13) (bottom). For a same system, $e_1$ and $e_2$ are plotted with the same color. Anti-Lagrange-like systems comply with $m_1e_1 = m_2e_2$ and both plots are easily distinguished while for Lagrange-like system, such that $e_1 = e_2$, both plots are almost overlaid. This is in agreement with what we deduce from Fig. 6.1 and Table 6.2.</figcaption></figure>

simulation of the complete Eqs. (C.13), respectively. We give these times in Table 6.3, for reference. In our model, the exterior of the co-orbital resonance does not exist (Fig. 3.1), and only the simulation of the complete Eqs. (C.13) allows $\tau_{\text{dest}}$, the time before destruction of the co-orbital pair, to be known. We also give this time in Table 6.3. The relative error between the three times $\tau_{\text{hs}}$ is consistently smaller than 1%, and the model is reliable on the whole tadpole region. For this choice of $\iota$, we have, for the first 5 systems, $\tau_{\text{hs}}/\tau_{\text{dest}} \approx 1.1$, which allowed us to predict the destruction time for system 6.

As expected from the values of $\tau_{\text{AL}}$ and $\tau_{\text{L}}$ given in Table 6.2, the three shortest simulations correspond to Lagrange-like systems, while the three longest correspond to anti-Lagrange like systems. This is particularly clear in Fig. 6.4, which displays the

<table><caption>Table 6.3 — Times to reach horseshoe-shaped orbits and until destruction, in millions of periods.</caption><thead><tr><th>#</th><th>τ<sub>hs</sub><sup>(a)</sup></th><th>τ<sub>hs</sub><sup>(s)</sup></th><th>τ<sub>hs</sub><sup>(c)</sup></th><th>|τ<sub>hs</sub><sup>(a)</sup> - τ<sub>hs</sub><sup>(s)</sup>|/τ<sub>hs</sub><sup>(s)</sup></th><th>|τ<sub>hs</sub><sup>(s)</sup> - τ<sub>hs</sub><sup>(c)</sup>|/τ<sub>hs</sub><sup>(c)</sup></th><th>τ<sub>dest</sub></th></tr></thead><tbody><tr><td>1</td><td>6.2753</td><td>6.2856</td><td>6.2427</td><td>0.001639</td><td>0.00683</td><td>6.9202</td></tr><tr><td>2</td><td>14.281</td><td>14.306</td><td>14.178</td><td>0.001748</td><td>0.00894</td><td>16.146</td></tr><tr><td>3</td><td>31.161</td><td>31.272</td><td>31.054</td><td>0.003550</td><td>0.00695</td><td>34.841</td></tr><tr><td>4</td><td>118.68</td><td>118.98</td><td>118.15</td><td>0.002521</td><td>0.00693</td><td>129.36</td></tr><tr><td>5</td><td>303.74</td><td>302.79</td><td>300.44</td><td>0.003137</td><td>0.00777</td><td>328.04</td></tr><tr><td>6</td><td>5467.9</td><td></td><td></td><td></td><td></td><td>~ 6.010</td></tr></tbody></table>
<figure><img src="image_19.png" /><figcaption>Fig. 6.6 — Total normalized energy of system 1 (**left**) and system 5 (**right**). The very anti-Lagrange-like system 5 features a pronounced plateau in its energy decrease and lives long while the very Lagrange-like system 1 does not show any marked plateau and has a short life.</figcaption></figure>

value of $\varpi_1 - \varpi_2$. For the short-lived systems (blue, green and red), the difference of the longitude of pericentres first settles around $60^\circ = \pi/3$ at low libration amplitude, before it moves to $0^\circ$, when in horseshoe orbit. For the long-lived systems (yellow, black and purple), the difference of the longitude of pericentres first settles around $240^\circ = 4\pi/3$ at low libration amplitude, before it moves to $180^\circ = \pi$, when in horseshoe orbit. This is in total agreement with Eqs. (3.44), (3.45), (3.49) & (3.50).

Figure 6.5 also confirms the Lagrange-like behaviour of short-lived systems, since $e_1 = e_2$ for these systems, while $m_1e_1 = m_2e_2$ for long-lived simulations, characteristic of their anti-Lagrange-like behaviour (Eqs. (3.44) & (3.45), Sect. 3.2.2). As also expected from the theoretical results, Lagrange-like system are still eccentric when they are old (e.g. the green and blue plots), while anti-Lagrange-like system are circular when they are at the end of their life (yellow plot).

### 6.3.3 A plateau in the energy decrease

In Fig. 6.6, we show the total energy $H_K + \iota H_P + T_1 + T_2$ of system 1 and system 5. The energy of the anti-Lagrange-like system 5 (yellow) features two steep decrease and a broad plateau. At the beginning of the simulation, the eccentricities of the planets are still significant, meaning that $\dot{w}_j$, the time derivative of their true longitude, is not constant. This prevents a solid rotation around the star, although the equilibrium rotation is already reached (small $\tau_{\text{rot}}^{(j)}$), and ensures dissipation: it is the first steep decrease. Then, the eccentricities are almost damped while the libration amplitude did not significantly increase yet. The motion around the star is almost a solid rotation and very few energy is dissipated into heat: it is the plateau. Finally, the high libration amplitude of $\xi$ reached in late tadpole and horseshoe-shaped orbits ensures the non-constancy of $\dot{w}_j$ and energy dissipation: it is the second steep decrease.
For almost restricted systems, very anti-Lagrange-like, the amplitude of libration has to be very high for the least massive planet to significantly perturb the solid rotation of the most massive one. Thus, these systems have very broad plateau and a very high life expectancy. On the other hand, the energy of the Lagrange-like system 1 (blue), does not feature such a marked plateau, because the eccentricities are always nonzero even at high libration amplitude, and hence its life is far shorter.

## 6.4 Tides raised on the central body

In the previous sections of this chapter, we discarded tides raised on the central body, because we showed that for a system star–planet–planet, the ratio of Eq. (6.1) is much smaller than unity. In this section, we show that such tides have to be taken into account for a system planet–satellite–satellite and we establish the corresponding tidal timescales.

### 6.4.1 Orders of magnitude

In Sect. 6.1.1, we only retained two out of the 12 tidal potentials of the form of Eq. (5.45), namely those corresponding to tides raised by the star on each planet and interacted with by the star. In this section, the central body is a planet, while the co-orbitals are satellites. Four additional contributions have hence to be taken into account, corresponding to tides raised by both satellites on the planet and interacted with by both satellites. Indeed, the ratio of Eq. (6.1) in that case can be of the order of unity.

The best example of such system is the pair Janus and Epimetheus orbiting Saturn. For this system, we have $l^{1/3} = 0.00163$ and according to Table 5.1, the $\kappa_2/Q$ of Saturn is $\kappa_2^{(0)}/Q_0 = 1.6 \times 10^{-4}$ (Lainey, 2016). The $\kappa_2/Q$ of Janus and Epimetheus are unknown, but we can estimate them. We first estimate the value of $\kappa_2$ with Eq. (5.31). For icy satellites made up of intact ice, it is common to take $\mu = 4$ GPa (Quillen et al., 2017), but fractures and pores reduce this value. According to Nimmo and Schenk (2006), the porous or fractured nature of the satellites is unlikely to decrease the value of $\mu$ by more than one order of magnitude (see their Fig. 6), so we consider $0.4 \text{ GPa} \le \mu \le 4 \text{ GPa}$ for both satellites and Eq. (5.31) yields$^{10}$ $1.12 \times 10^{-5} \le \kappa_2 \le 1.12 \times 10^{-4}$ for Epimetheus and $2.74 \times 10^{-5} \le \kappa_2 \le 2.74 \times 10^{-4}$ for Janus. These values are rough approximations and we consider, for both satellites, $\kappa_2 = 10^{-4}$ if $\mu = 0.4$ GPa and $\kappa_2 = 10^{-5}$ if $\mu = 4$ GPa.

We now give a minimal value for the quality factor $Q$. According to the considerations stated in the last paragraph of Sect. 5.2.3, the time-lag $\Delta t$ cannot exceed the time needed for a seismic wave to travel through the satellite$^{11}$. Calling $K$ the bulk modulus, the speed of sound is given by $c_s = ((K + 4\mu/3)/\rho)^{1/2}$. Assuming $K \sim \mu$, the time lag $\Delta t$ is majored by $\Delta t \le R_j/c_s = R_j(3\rho/(7\mu))^{1/2}$ and the quality factor $Q$ complies with (Eq. (5.44))

$$Q \ge \eta^{-1} R_j^{-1} \sqrt{\frac{7\mu}{3\rho}}. \qquad (6.42)$$

<footnote>$^{10}$Using $\rho = 618 \text{ kg/m}^3$, $R_j = 58.1 \text{ km}$ and $g = 0.0104 \text{ m/s}^2$ for Epimetheus and $\rho = 646 \text{ kg/m}^3$, $R_j = 89.5 \text{ km}$ and $g = 0.0158 \text{ m/s}^2$ for Janus.</footnote>

<footnote>$^{11}$In practice, $\Delta t$ is much smaller than that.</footnote>
Assuming a typical frequency of excitation $\eta \sim 2$ rad/day, $\mu = 0.4$ GPa gives $Q \ge 914$ for Epimetheus and $Q \ge 580$ for Janus, while $\mu = 4$ GPa yields $Q \ge 2890$ for Epimetheus and $Q \ge 1835$ for Janus. These are once again very rough estimates and we consider for both satellites $Q \ge 1000$ if $\mu = 0.4$ GPa and $Q \ge 3000$ if $\mu = 4$ GPa. We thus have for Janus and Epimetheus, and more generally for any satellite far from the hydrostatic equilibrium

$$ \frac{\kappa_2^{(j)}}{Q_j} \le \begin{cases} 10^{-7} & \text{if } \mu = 0.4 \text{ GPa,} \\ 0.3 \times 10^{-8} & \text{if } \mu = 4 \text{ GPa.} \end{cases} \qquad (6.43) $$

In the case of the system Saturn–Janus–Epimetheus, the ratio of Eq. (6.1) verifies

$$ l^{1/3} \frac{\kappa_2^{(0)}}{Q_0} \ge \begin{cases} 2.6 & \text{if } \mu = 0.4 \text{ GPa,} \\ 87 & \text{if } \mu = 4 \text{ GPa.} \end{cases} \qquad (6.44) $$

and tides raised on the central planet cannot be discarded since this ratio is not much smaller than 1. It is unclear if it is much larger than 1 for this system, so we cannot say if the two tidal contributions that we considered in the previous sections of this chapter (tides raised on the orbiting bodies by the central body and interacted with by the central body) can be discarded in this section. If they cannot, then the real parts given by Eqs. (6.28) have to be added to the real parts given by Eqs. (6.61) in order to get to total tidal contribution.

For a system composed of two large satellites (like the Galilean satellites) orbiting a gas giant, $\kappa_2^{(j)}/Q_j$ is of the order of 0.015 (see the value for Io in Table 5.1), while $\kappa_2^{(0)}/Q_0$ is rather of the order of $10^{-5}$ to $10^{-4}$ (see the values for Jupiter and Saturn in Table 5.1). The ratio (6.1) is hence much smaller than 1 and tides raised on the central planet can be discarded again, as if it was a star–planet–planet system, meaning that the results of Sects. 6.2 & 6.3 apply.

### 6.4.2 Expression of the pseudo-Hamiltonian

As in Sect. 6.1.1, we write the tidal perturbation of the pseudo-Hamiltonian as

$$ H_t = U_t^{(11)} + U_t^{(22)} + U_t^{(12)} + U_t^{(21)} + T, \qquad (6.45) $$

where $U_t^{(jk)}$ is the contribution of tides raised by the satellite $k$ on the planet and interacted with by the satellite $j$, and $T$ is the kinetic energy of rotation of the planet. We denote $R_0$ the radius of the planet and its rheology is characterized by its second Love number $\kappa_2^{(0)}$ (Eqs. (5.35) & (5.36)) and its quality factor $Q_0$ (Sect. 5.2.4). We have (Eq. (5.45))

$$ U_t^{(jk)} = -\kappa_2^{(0)} G m_j m_k \frac{R_0^5}{r_j^3 r_k^{\star 3}} P_2(\cos S), \quad T = \frac{\Theta'^2}{2C_0} \qquad (6.46) $$

with

$$ S = \lambda_j - \lambda_k - (\theta - \theta^*), \qquad (6.47) $$

where $\theta$ is the rotation angle of the planet, $\Theta' = C_0\omega$ is the conjugated momentum of $\theta$, $C_0 = \alpha_0 m_0 R_0^2$ is the moment of inertia of the planet, $\omega = d\theta/dt$ is its rotation rate, and
$\alpha_0$ is a dimensionless structure constant depending on the state equation of the planet ($\alpha_j = 2/5$ if it is homogeneous). Like in Sect. 6.1.1, we denote

$$z^* = z(t - \Delta t_0) \qquad (6.48)$$

for any quantity $z$. We repeat the calculations from Sect. 6.1.1 and we write the tidal Hamiltonian $H_t$ in the variables $(J, J_2, \xi, \xi_2, X_j, \bar{X}_j)$. We define $\Theta = \Theta'/(m\bar{a}^2\eta)$ in order to normalise $\Theta'$ (see Eq. (3.18)). We average $U_t^{(jk)}$ over $\lambda_j$ and $\lambda_k$ and $U_t^{(jk)}$ depends on these angles through $\Delta\lambda^{(jk)}$ only, where

$$\Delta\lambda^{(jk)} = \lambda_j - \lambda_k^* \quad \text{and} \quad \Delta\theta = \theta - \theta^*. \qquad (6.49)$$

In virtue of Eq. (6.47), $U_t^{(jk)}$ depends on $\theta$ through $\Delta\theta$ only. We define the quantities

$$q_0 = \kappa_2^{(0)} \varphi_0^5, \quad \text{and} \quad \varphi_0 = \frac{R_0}{\bar{a}}, \qquad (6.50)$$

and performing the rescaling (3.22), we get for $\mathcal{H}_t = H_t/(m\bar{a}^2\eta)$ the expression<fnref n="12" />

$$\begin{gathered} \mathcal{U}_t^{(jk)} = -q_0\eta \frac{m_j m_k}{m_0 m} \mathcal{R}_j^{-6} \mathcal{R}_k^{*-6} \left\{ A_t^{(jk)} + \mathcal{D}_2^{(jk)} + \mathcal{D}_4^{(jk)} \right\} = \frac{U_t^{(jk)}}{m\bar{a}^2\eta}, \qquad (6.51) \\ \mathcal{T} = \frac{\eta}{2\alpha_0} \frac{m}{m_0} \frac{\Theta^2}{\varphi_0^2} = \frac{T}{m\bar{a}^2\eta}, \end{gathered}$$

where (Eq. (6.9))

$$\mathcal{D}_2^{(jk)} = B_t^{(jk)} (\mathcal{R}_j^{-1} X_j \bar{X}_j + \mathcal{R}_k^{*-1} X_k^* \bar{X}_k^*) + (\mathcal{R}_j \mathcal{R}_k^*)^{-1/2} (C_t^{(jk)} X_j \bar{X}_k^* + \bar{C}_t^{(jk)} X_k^* \bar{X}_j), \quad (6.52)$$

and the coefficients $A_t^{(jk)}$, $B_t^{(jk)}$ and $C_t^{(jk)}$ are given by Eq. (6.10), substituting $\Delta\lambda^{(jk)}$ for $\Delta\lambda_j$ and $\Delta\theta$ for $\Delta\theta_j$. The coefficient $\mathcal{D}_4^{(jk)}$ contains the fourth order in eccentricity and is given in appendix B.3, with the same substitution. As in Sect. 6.1.1, the pseudo-Hamiltonian that we consider for the model is

$$\mathcal{H} = \mathcal{H}_K(J, J_2) + i\mathcal{H}_P(\xi, X_1, X_2, \bar{X}_1, \bar{X}_2) + \mathcal{H}_t, \qquad (6.53)$$

where $\mathcal{H}_K$ is given by Eq. (6.12), the perturbation $i\mathcal{H}_P = \mathcal{H}^{(0)} + \mathcal{H}^{(2)} + \mathcal{H}^{(4)}$ is given by Eq. (3.29) and

$$\mathcal{H}_t = \sum_{j=1}^{2} \sum_{k=1}^{2} \mathcal{U}_t^{(jk)} (J, J_2, J^*, J_2^*, \xi, \Delta\xi, \Delta\xi_2, \Delta\theta, \mathbf{X}, \mathbf{X}^*, \bar{\mathbf{X}}, \bar{\mathbf{X}}^*) + \mathcal{T}(\Theta), \qquad (6.54)$$

where $\mathbf{X} = (X_1, X_2)$, $\Delta\xi = \xi - \xi^*$ and $\Delta\xi_2 = \xi_2 - \xi_2^*$.

### 6.4.3 Equilibria and eigenvalues

As in Sect. 6.1.2, we derive the equations of motion from the pseudo-Hamiltonian (6.53), following the procedure described in Sect. 5.3 and using Eq. (3.23). The major difference with Sect. 6.1 is the existence of crossed tides. While the Hamiltonian $\mathcal{H}_t$ depends only

<footnote n="12">¹²$\mathcal{R}_j(J, J_2)$ is defined by Eq. (3.20) and $m = \sqrt{m_1 m_2}$.</footnote>
on $\lambda_j - \lambda_j^*$ in Sect. 6.1, we now also have a dependency on $\lambda_j - \lambda_k^*$, where $j \neq k$. As a consequence, $H_t$ depends on $\xi$ in Eq. (6.54), whereas it did not depend on $\xi$ in Eq. (6.13). This new dependency complexifies immensely the differential system derived from the pseudo-Hamiltonian, and while we were able to give it in a compact form in Eq. (6.15), it is no longer possible to write the differential system in a readable form. Instead, we merely linearize it in the vicinity of its equilibria.

Due to the conservation of the total angular momentum, the differential system derived from the pseudo-Hamiltonian (6.53) only provides five independent equations for six variables to locate the equilibria, and in order to have one unique equilibrium in the vicinity of which we linearize the differential system, we again choose the arbitrary Eq. (6.22). As explained in Sect. 6.2.1, the linearized system does not depend on the choice of Eq. (6.22) at first order in $\iota q_0$. We define $\vartheta = 1 - \omega/\eta$ and the Lagrangian equilibrium $L_4$ is located at<fnref n="13" />

$$ 
\begin{aligned}
\vartheta_0 &= X_{1,0} = X_{2,0} = 0, \\
\xi_0 &= \frac{\pi}{3} + \frac{q_0}{\sqrt{3}}, \\
a_{1,0} &= \bar{a} \left( 1 + 2q_0 \frac{m_1 - m_2}{m_0} \right), \\
a_{2,0} &= \bar{a} \left( 1 + 2q_0 \frac{m_2 - m_1}{m_0} \right).
\end{aligned} 
\qquad (6.55) 
$$

Let $\mathcal{X}_0 = {}^t(J_{2,0}, \vartheta_0, J_0, \xi_0, X_{1,0}, X_{2,0})$ and $\mathcal{X} = {}^t(J_2, \vartheta, J, \xi, X_1, X_2)$. Linearized in the vicinity of the Lagrangian equilibrium $L_4$, the differential system derived from the pseudo-Hamiltonian (6.53) reads

$$ \frac{d}{dt} \Delta \mathcal{X} = (\mathcal{Q}'_0 + \mathcal{Q}'_1) \Delta \mathcal{X}, \qquad (6.56) $$

where $\Delta \mathcal{X} = \mathcal{X} - \mathcal{X}_0$. The matrix $\mathcal{Q}'_0$ derives from the conservative Hamiltonian $\mathcal{H}_K + i\mathcal{H}_P$, while $\mathcal{Q}'_1$ corresponds to the tidal contribution $\mathcal{H}_t$. Because the dynamics of $X_1$ and $X_2$ is uncoupled from the dynamics of the other variables, the matrices $\mathcal{Q}'_0$ and $\mathcal{Q}'_1$ are block diagonal and take the form

$$ \mathcal{Q}'_0 = \begin{pmatrix} \mathcal{Z}'_0 & 0_{4 \times 2} \\ 0_{2 \times 4} & \mathcal{M}'_0 \end{pmatrix} \quad \text{and} \quad \mathcal{Q}'_1 = \begin{pmatrix} \mathcal{Z}'_1 & 0_{4 \times 2} \\ 0_{2 \times 4} & \mathcal{M}'_1 \end{pmatrix}, \qquad (6.57) $$

where $\mathcal{Z}'_0, \mathcal{M}'_0, \mathcal{Z}'_1$ and $\mathcal{M}'_1$ are given in this Maxima worksheet. The set of eigenvalues of $\mathcal{Q}_0$ is (Eq. (6.25))

$$ \{0, 0, i\nu, -i\nu, ig_1, ig_2\}, \qquad (6.58) $$

where $\nu$ (Eq. (3.32)) is the libration frequency of $\xi$ around $L_4$ while $g_1$ and $g_2$ (Eq. (3.42)) are the precession frequencies of the pericentres in the eccentric eigenmodes anti-Lagrange and Lagrange, respectively. The first 0 eigenvalue corresponds to the constancy of $J_2$ (conservation of the total angular momentum), while the second one corresponds to the constant rotation rate of the planet. Like in Sect. 6.2.1, all six eigenvalues of $\mathcal{Q}'_0$ are pure imaginary and the linearized system is quasiperiodic without tides.

<footnote n="13">¹³Or equivalently, $J_0 = \frac{9}{4}q_0 \frac{m(m_1 - m_2)}{m_0(m_1 + m_2)}$ and $J_{2,0} = \frac{9}{8}q_0 \frac{(m_1 - m_2)^2}{m_0m}$.</footnote>
We use the perturbative approach described in Sect. 2.4.2 to find the eigenvalues of $Q'_0 + Q'_1$ at first order in $\iota q_0$. In order to alleviate the expression of the eigenvalues, it is convenient to define, for integers $p$ and $q$, the quantity

$$ \mathfrak{N}_p^q = \frac{p(m_1 - m_2)^2 + qm_1m_2}{m_0(m_1 + m_2)}. \qquad (6.59) $$

The set of eigenvalues of $Q'_0 + Q'_1$ is

$$ \{\lambda, 0, \lambda, \bar{\lambda}, \lambda_{\text{AL}}, \lambda_{\text{L}}\}, \qquad (6.60) $$

with

$$ \begin{aligned} \lambda &= \boxed{3\eta \frac{q_0}{Q_0} \mathfrak{N}_1^1 (3 - \alpha_0^{-1} \vartheta_0^{-2} \iota)}, \\ \lambda &= \boxed{\frac{27}{2} \eta \frac{q_0}{Q_0} \mathfrak{N}_0^1} + i\nu - i\nu \frac{q_0}{8} \mathfrak{N}_7^{-206} - \frac{45}{16} i\eta^2 \nu^{-1} q_0 \iota, \\ \lambda_{\text{AL}} &= \boxed{-\frac{225}{8} \eta \frac{q_0}{Q_0} \mathfrak{N}_0^1} + ig_1 + \frac{3}{8} i\eta q_0 \mathfrak{N}_7^{45}, \\ \lambda_{\text{L}} &= \boxed{-\frac{3}{8} \eta \frac{q_0}{Q_0} \mathfrak{N}_{28}^{37}} + ig_2 + \frac{15}{8} i\eta q_0 \mathfrak{N}_4^7. \end{aligned} \qquad (6.61) $$

The eigenvalues are no longer pure imaginary and we boxed the real parts. They are proportional to the inverse quality factor $Q_0^{-1}$ and as in Sect. 6.2.1, elastic tides do not contribute to the real parts but only change slightly the fundamental frequencies. Due to the conservation of the total angular momentum, the matrix $Q'_0 + Q'_1$ is singular and one of its eigenvalues is zero. The eigenvalues $\lambda_{\text{AL}}$ and $\lambda_{\text{L}}$, perturbations of $ig_1$ and $ig_2$, respectively, have non-zero negative real parts and both eccentric eigenmodes Lagrange and anti-Lagrange are damped to zero, as with the tides raised on the co-orbitals. Similarly, the real part of $\lambda$ and $\bar{\lambda}$, perturbations of $i\nu$ and $-i\nu$, are strictly positive, which once again leads to an exponential increase of the libration amplitude of $\xi$ when the system is close to the Lagrangian equilibria $L_4$ or $L_5$.

A notable difference with Sect. 6.2.1 resides in the eigenvalue $\lambda$, associated with the rotation of the central planet. While $\lambda_1$ and $\lambda_2$ were large negative in Sect. 6.2.1 (Eq. (6.28)), leading to a quick damping of the rotation rates of the co-orbitals, $\lambda$ is not necessarily negative, and can be of the same order of magnitude as the other eigenvalues' real parts. As a consequence, the rotation of the central planet does not necessarily evolve on timescales much smaller than the timescale of eccentricity damping and libration amplitude pumping. If $\lambda < 0$, then the system converges towards the spin-orbit synchronization<fnref n="14" />, while $\lambda > 0$ means that the system diverges away from the spin-orbit synchronization. There exists a critical semimajor axis $\bar{a}_c$ beyond which (resp. below which) $\lambda < 0$ and the system is stable (resp. $\lambda > 0$ and the system is unstable). According to the value of $\lambda$ in Eq. (6.61), such critical semimajor axis is given by

$$ \bar{a}_c = \sqrt{3\alpha_0\iota^{-1}R_0}. \qquad (6.62) $$

<footnote n="14">¹⁴Or towards the spin-orbit pseudo-synchronization at non-zero eccentricity.</footnote>
This result was reached in a completely different manner by Hut (1980). Hut studied a binary system and showed that the Hessian of the Hamiltonian at the equilibrium (circularity, coplanarity and spin-orbit synchronization) is positive-definite if, and only if, $\bar{a} \ge \bar{a}_c$, hence the result<fnref n="15" />. Hut's result can be reformulated equivalently by saying that the spin-orbit synchronization is stable if, and only if, the central planet contains in its rotation at most 1/4 of the system's total angular momentum. In the case of tides raised on the co-orbitals, the eigenvalues $\lambda_1$ and $\lambda_2$ in Eq. (6.28) are always negative because the co-orbitals, due to their low masses and radii, can never contain in their rotations more than 1/4 of the total angular momentum.

<figure><img src="image_20.png" /><figcaption>Fig. 6.7 — Schematic representation of the trajectories of the co-orbital satellites in the case of a non-synchronized planet. We distinguish the cases $\bar{a}_s < \bar{a}_c$ (left) and $\bar{a}_s > \bar{a}_c$ (right). For visibility, the spiraling motion is largely exaggerated and only one co-orbital satellite is represented.</figcaption></figure>

When the ratio $\omega/\eta$ is below 1, then the co-orbitals are in advance (by an angle $(\eta - \omega)\Delta t_0$) with respect to the tidal bulges they raise on the central planet. The subsequent torque increases $\omega$ and, by conservation of the total angular momentum, decreases the semimajor axes. When $\omega/\eta > 1$, $\omega$ decreases while the semimajor axes increase. The ratio $\omega/\eta$ is below 1 (resp. above 1) if $\bar{a}$ is smaller (resp. greater) than the synchronous semimajor axis $\bar{a}_s$, given by

$$ \bar{a}_s = \left( \frac{Gm_0}{\omega^2} \right)^{1/3} . \qquad (6.63) $$

Four cases are possible :

* $\bar{a} \le \bar{a}_c$ and $\bar{a} \le \bar{a}_s$. The planet's rotation rate is sub-synchronous and $\lambda > 0$: The co-orbitals fall onto the planet.

<footnote n="15">¹⁵In Hut's work, the expression of $\bar{a}_c$ features $m_0/m_1$ instead of $\iota^{-1} = m_0/(m_1 + m_2)$.</footnote>
* $\bar{a} \le \bar{a}_c$ and $\bar{a} \ge \bar{a}_s$. The planet's rotation rate is super-synchronous and $\lambda > 0$: The co-orbitals go to infinity.

* $\bar{a} \ge \bar{a}_c$ and $\bar{a} \le \bar{a}_s$. The planet's rotation rate is sub-synchronous and $\lambda < 0$: The co-orbitals' semimajor axes decrease and stabilize.

* $\bar{a} \ge \bar{a}_c$ and $\bar{a} \ge \bar{a}_s$. The planet's rotation rate is super-synchronous and $\lambda < 0$: The co-orbitals' semimajor axes increase and stabilize.

The situation is illustrated in Fig. 6.7. For the system Saturn–Janus–Epimetheus, we consider $\alpha_0 = 0.22$ (Fortney et al., 2018, Sect. 3.4.2) and a numerical evaluation yields $\bar{a}_c = 1.2 \times 10^4 R_0$. The critical semimajor axis $\bar{a}_c$ is larger than 12000 Saturn's radii and Janus and Epimetheus are largely below (The rotation of Saturn contains 99.999992 % of the total angular momentum of the system Saturn–Janus–Epimetheus). Since Saturn is super-synchronous ($\omega/\eta = 1.5813$), we expect Janus and Epimetheus' semimajor axes to diverge to infinity under the influence of the tides they raise on Saturn (green plot in Fig. 6.7). Deimos around Mars is also in the green case while Phobos is in the black case and the Moon around the Earth is in the orange case<fnref n="16" />.

### 6.4.4 Characteristic timescales

If, for the considered co-orbital system, the ratio (6.44) is found much larger than unity, then the eigenvalues are given by Eq. (6.61) only. If it is of the order of unity, the eigenvalues given by Eq. (6.28) have to be added to the eigenvalues given by Eq. (6.61). Finally, if the ratio (6.44) is much smaller than 1, only the contributions of Eq. (6.28) need to be considered and the timescales are given by Eq. (6.29). Once the eigenvalues are established, the timescales, defined as the times needed for the associated eigenmodes to be damped or pumped by a factor $\exp(1)$, are given by

$$ \tau_{\text{lib}} = \frac{T}{2\pi |\mathfrak{Re} \, \lambda|}, \quad \tau_{\text{AL}} = \frac{T}{2\pi |\mathfrak{Re} \, \lambda_{\text{AL}}|}, \quad \tau_{\text{L}} = \frac{T}{2\pi |\mathfrak{Re} \, \lambda_{\text{L}}|}, \quad \tau_{\text{rot}} = \frac{T}{2\pi |\mathfrak{Re} \, \lambda|}. \quad (6.64) $$

We have $\tau_{\text{lib}} = 25\tau_{\text{AL}}/12$ and as in Sect. 6.2.2, $\tau_{\text{lib}}$ is always of the same order as $\tau_{\text{AL}}$ and anti-Lagrange-like system are nearly circular when old. While the ratio $\tau_{\text{AL}}/\tau_{\text{L}}$ depends on $x = m_1/m_2$ and $y = q_2Q_1/(q_1Q_2)$ in Sect. 6.2.2, where tides were raised on the co-orbitals, it now depends on $x$ alone and reads

$$ \frac{\tau_{\text{AL}}}{\tau_{\text{L}}} = \frac{37}{75} + \frac{28}{75} \frac{(1-x)^2}{x}. \quad (6.65) $$

In Fig. 6.8, we plot the ratio $\tau_{\text{AL}}/\tau_{\text{L}}$. The equation $\tau_{\text{AL}} = \tau_{\text{L}}$ has roots $\{x_1, x_2\} = \{(47 \pm 5\sqrt{57})/28\} \approx \{0.3304, 3.0268\}$. Between the roots, $\tau_{\text{AL}} \le \tau_{\text{L}}$ and the system is Lagrange-like, while it is anti-Lagrange-like outside the roots (Sect. 3.2.2). In Sect. 6.2.2, there exists values of $x$ and $y$ such that $\tau_{\text{AL}} \ll \tau_{\text{L}}$ (red-black regions in Fig. 6.1), and the associated systems are very Lagrange-like and short-lived. Here, the ratio $\tau_{\text{AL}}/\tau_{\text{L}}$ is never smaller than $37/75 \approx 1/2$ and co-orbital systems are never very Lagrange-like when tides are raised on the central body.

<footnote n="16">¹⁶Because of the Sun's perturbations, the Moon will still leave the Earth's Hill sphere.</footnote>
<figure><img src="image_21.png" /><figcaption>Fig. 6.8 — Value of the ratio $\tau_{AL}/\tau_{L}$ as a function of $x = m_{1}/m_{2}$. Between (resp. outside) the roots of $\tau_{AL}/\tau_{L} = 1$, the system settles in the Lagrange (resp. anti-Lagrange) configuration.</figcaption></figure>

Given the value of $\lambda$ in Eq. (6.61) from Sect. 6.4.3, the timescale $\tau_{\text{rot}}$ of convergence to, or divergence away from, the spin-orbit synchronization is

* Of the same order of magnitude as the other tidal timescales if $\bar{a} \ll \bar{a}_c$.

* Much larger than the other tidal timescales if $\bar{a} \approx \bar{a}_c$.

* Much smaller than the other tidal timescales if $\bar{a} \gg \bar{a}_c$.

## 6.5 Conclusion

In this chapter, we extended the work of Chap. 3 by considering tidal dissipation in our model. We first proved that for a co-orbital system star–planet–planet, it is enough to only consider tides raised by the star on the planets and interacted with by the star. We built the pseudo-Hamiltonian of the co-orbital dynamics, including these types of tides, and we gave the complete expression of the averaged equations of motion, at fourth order in eccentricity. We provided three independent verifications of the differential system (6.15), namely, the conservation of the total angular momentum (Sect. 6.1.3), the pseudo-synchronization (Eq. 6.23) and the numerical simulations (Sect. 6.3).

We analytically computed the eigenvalues of the system in the vicinity of the Lagrangian equilibria and we showed that these equilibria are made unstable by tides, leading to the systematic destruction of the system by close encounters. The tidal evolution is essentially dominated by three characteristics timescales. Two of them are responsible for damping the proper modes of the eccentricity and the third one for pumping the libration amplitude of the resonant angle $\lambda_1 - \lambda_2$. The lifetime of the system depends on the
masses, radii, tidal parameters and semimajor axes, and we provided a precise analytical expression that we applied to hypothetical cases in Table 6.1 and to the detection of co-orbital exoplanets in Sect. 6.2.3. From the computation of the eigenvalues, we also concluded that, depending on the parameters, the system can settle along the Lagrange configuration, or along the anti-Lagrange one. The latter configuration is associated with a long lifetime, while the former configuration is associated with a short lifetime.

The validity of the analytical work was established in Sect. 6.3, where we performed numerical simulations of targeted systems, using our model as well as the complete equations of motion of the *N*-body problem. These simulations show the great reliability of the model, and due to the average over the mean motions, it is $\eta/\nu$ times faster to integrate than the complete equations with similar CPU performances (27 times faster with $\iota = 2 \times 10^{-4}$), while also generating less oscillating curves.

The main achievement of this chapter was to provide an analytical proof of the instability of the Lagrangian equilibria under tides and to give the timescale of destruction. Recently, Dobrovolskis and Lissauer (2022) gave another proof with another model. They consider a restricted problem where one of the co-orbitals, the only one where tides are raised, is massless and point-mass. Although this approach does not allow them to provide a destruction timescale and is limited to the restricted case, it allows them to prove the instability of the Lagrangian equilibria without average over the mean longitudes nor expansion in series of the eccentricities or around the Keplerian resonance.

Finally, we dedicated Sect. 6.4 to the case where tides raised by the co-orbitals on the central body cannot be discarded, as it is the case with a system planet–satellite–satellite made up of small satellites. We provided a complete discussion that allowed us to determine which types of tides need to be taken into account, depending on the type of system. We computed the eigenvalues in the vicinity of the Lagrangian equilibria in the case of tides raised on the central body, and we inferred the associated timescales of tidal evolution. We showed that these tides have a lot of similarities with those raised on the orbiting bodies, as well as some proper characteristics, like the existence of a critical semimajor axis that decides whether or not the spin-orbit synchronization is stable. This last section is important to understand the stability of the co-orbital satellites of Saturn and it will be explored in more details in a forthcoming study.
# Chapter 7

## Tides in the planar four-body
$p:p:p+1$ resonance chain

*Most results of this chapter were first published in Couturier et al. (2022).*

### 7.1 Pseudo-Hamiltonian and equations of motion

In this chapter, we push further the work of Chap. 4 on the planar resonance chain $p:p:p+1$ by removing the point-mass approximation. Tidal dissipation due to the distortion of the extended bodies is hence responsible for long-term changes on the dynamics of the chain. We want here to investigate the effects of tides following the same procedure that we used in Chap. 6, where we studied tides in the resonance chain 1:1. More precisely, we want to

* Establish the pseudo-Hamiltonian of the averaged problem,

* Derive the equations of motion from the pseudo-Hamiltonian,

* Locate the equilibria of the averaged problem and linearize in their vicinity,

* Obtain stability and timescales results from the eigenvalues,

* Check the results given by the eigenvalues with numerical simulations.

The star and three planets have radii $R_0$ to $R_3$, respectively, and we use the same notations as in Chaps. 4 & 6.

#### 7.1.1 Pseudo-Hamiltonian

We study a four-body system star-planet-planet-planet and as we justified in Sect. 6.1.1, among the $4(4-1)^2 = 36$ tidal contributions as given by Eq. (5.45), only three
of them need to be retained. We thus write the perturbation to the Hamiltonian due to tides as (Eq. (6.2))

$$H_t = \sum_{j=1}^{3} U_t^{(j)} + \sum_{j=1}^{3} T_j, \qquad (7.1)$$

where the tidal potential $U_t^{(j)}$ and the kinetic energy of rotation $T_j$ of planet $j$ are given by Eq. (6.3), conserving the same notations<fnref n="1" />. In particular, we reuse the notations $z^* = z(t - \Delta t_j)$ and $\Delta\zeta = \zeta - \zeta^*$ defined by Eqs. (6.5) & (6.6). We also denote (Eq. (6.7)) $q_j = \kappa_2^j \varphi_j^5$ and in this chapter we have

$$\varphi_j = \frac{R_j}{a_j^{(0)}}, \qquad (7.2)$$

which coincides with the definition of Chap. 6, $\varphi_j = R_j/\bar{a}$, for $j = 1, 2$. We define $\mathcal{R}_j$ similarly as in Eq. (3.20), that is

$$\mathcal{R}_j = \left( \frac{a_j}{a_j^{(0)}} \right)^{1/2} = \frac{\tilde{\Lambda}_j}{\Lambda_j^*} = C_j \Lambda_j, \qquad (7.3)$$

where $C_j$ is defined by Eq. (4.14), $\Lambda_j^*$ by<fnref n="2" /> Eq. (4.5) and $\Lambda_j$ by Eq. (4.10). In Sect. 6.1.1, we wrote the tidal Hamiltonian $H_t$ in the variables $(J, J_2, X_j; \xi, \xi_2, \bar{X}_j)$, but to be consistent with the work of Chap. 4, we write it here in the variables $(D_j, L, \Gamma, G; \sigma_j, \xi, \xi_2, \xi_3)$, introduced in Sect. 4.1.1. As in Sect. 6.1.1, we normalize the action variable $\Theta_j'$ by writing $\Theta_j = \Theta_j'/\Gamma^*$, where the constant $\Gamma^*$ is given by Eq. (4.9). In order to only consider the secular (i.e. long-term) dynamics, we average $U_t^{(j)}$ over $\lambda_j$ and $U_t^{(j)}$ depends on this angle through $\Delta\lambda_j$ only. Similarly, $U_t^{(j)}$ depends on $\theta_j$ through $\Delta\theta_j$ only. We perform the rescaling given by Eq. (4.10) and at second order in eccentricity, the tidal Hamiltonian $\mathcal{H}_t = H_t/\Gamma^*$ reads

$$\begin{align}
\mathcal{U}_t^{(j)} &= -q_j \frac{m_0}{m_j} n_j^{(0)} C_j^{-1} \mathcal{R}_j^{-6} \mathcal{R}_j^{*-6} (A_t^{(j)} + \mathcal{D}_2^{(j)}) = \frac{U_t^{(j)}}{\Gamma^*}, \tag{7.4} \\
\mathcal{T}_j &= \frac{C_j n_j^{(0)} \Theta_j^2}{2 \alpha_j \varphi_j^2} = \frac{T_j}{\Gamma^*}, \nonumber
\end{align}$$

where (Eq. (6.9))

$$\mathcal{D}_2^{(j)} = B_t^{(j)} (\mathcal{R}_j^{-1} D_j + \mathcal{R}_j^{*-1} D_j^*) + (\mathcal{R}_j \mathcal{R}_j^*)^{-1/2} \mathfrak{Re} (C_t^{(j)} \sqrt{D_j D_j^*} e^{i \Delta \varpi_j}), \qquad (7.5)$$

and the coefficient $A_t^{(j)}$, $B_t^{(j)}$ and $C_t^{(j)}$ are given by Eq. (6.10). The total pseudo-Hamiltonian that we consider for the model is then (Eq. (6.11))

$$\mathcal{H} = \mathcal{H}_K(D_j, L, \Gamma, G) + i\mathcal{H}_P(D_j, \sigma_j, \xi) + \mathcal{H}_t, \qquad (7.6)$$

<footnote n="1">¹Do not confuse $C_j$, defined below Eq. (6.4), with $C_j$, defined by Eq. (4.14).</footnote>

<footnote n="2">²For the variables $\Lambda_j^*$ and $\Gamma^*$, the superscripted star does not denote the evaluation at time $t - \Delta t_j$ but the value at the Keplerian resonance given by Eq. (4.5).</footnote>
where $\mathcal{H}_K$ is given by Eq. (4.13), the perturbation $\iota \mathcal{H}_P = \mathcal{H}^{(0)} + \mathcal{H}_{1,3}^{(1)} + \mathcal{H}_{2,3}^{(1)} + \mathcal{H}_{1,2}^{(2)} + \mathcal{H}_{1,3}^{(2)} + \mathcal{H}_{2,3}^{(2)}$ by Eqs. (4.18), (4.19) & (4.20) and (Eq. (6.13))

$$ \mathcal{H}_t = \sum_{j=1}^{3} \mathcal{U}_t^{(j)}(\mathcal{R}_j, \mathcal{R}_j^*, \Delta\xi, \Delta\xi_2, \Delta\xi_3, \Delta\sigma_j, \Delta\theta_j) + \sum_{j=1}^{3} \mathcal{T}_j(\Theta_j). \quad (7.7) $$

The quantity $\mathcal{R}_j$ depends on $\Lambda_j$, that is on $L$, $\Gamma$, $G$ and the $D_j$.

### 7.1.2 Equations of motion

The equations of motion are derived from the pseudo-Hamiltonian (7.6) using the Hamilton Eqs. (2.18), following the procedure described in Sect. 5.3. The perturbation to the vector field $F_0$ (derived from the conservative Hamiltonian of Eq. (4.21)), due to tides and at second order in eccentricity, reads

$$
\begin{align}
\dot{D}_j &= 3n_j^{(0)} D_j \frac{m_0}{m_j} \mathcal{R}_j^{-13} \frac{q_j}{Q_j} (12\omega_j + 57\mathcal{R}_j - 76), \nonumber \\
\dot{L} &= 3\eta \frac{m_0}{m_1} \mathcal{R}_1^{-13} \frac{q_1}{Q_1} \{\omega_1 (\Lambda_1 + 27D_1) + (3\mathcal{R}_1 - 4) (\Lambda_1 + 46D_1)\}, \nonumber \\
\dot{\sigma}_j &= -\frac{15}{2} n_j^{(0)} \frac{m_0}{m_j} q_j \mathcal{R}_j^{-13}, \nonumber \\
\dot{\xi} &= \frac{3}{2} \eta \frac{m_0}{m_1} q_1 \mathcal{R}_1^{-13} \frac{4\Lambda_1 + 65D_1}{\Lambda_1} - \frac{3}{2} \eta \frac{m_0}{m_2} q_2 \mathcal{R}_2^{-13} \frac{4\Lambda_2 + 65D_2}{\Lambda_2}, \nonumber \\
\dot{\Gamma} &= \sum_{j=1}^3 \frac{3 n_j^{(0)}}{n_3^{(0)}} n_j^{(0)} \frac{m_0}{m_j} \mathcal{R}_j^{-13} \frac{q_j}{Q_j} \{\omega_j (\Lambda_j + 27D_j) + (3\mathcal{R}_j - 4) (\Lambda_j + 46D_j)\}, \nonumber \\
\dot{G} &= \sum_{j=1}^3 \frac{3 n_j^{(0)}}{m_j} \mathcal{R}_j^{-13} \frac{q_j}{Q_j} \{\omega_j (\Lambda_j + 15D_j) + (3\mathcal{R}_j - 4) (\Lambda_j + 27D_j)\}, \nonumber \\
\dot{\omega}_j &= -3n_j^{(0)} C_j \alpha_j^{-1} \gamma_j^{-2} \frac{m_0}{m_j} \mathcal{R}_j^{-13} \frac{q_j}{Q_j} \{\omega_j (\Lambda_j + 15D_j) + (3\mathcal{R}_j - 4) (\Lambda_j + 27D_j)\}, \nonumber
\end{align}
$$

where we posed $\omega_j = \omega'_j/n_j^{(0)}$. Since the differential system does not depend on $\xi_2$ and $\xi_3$ and as the dynamics of these fast-circulating angles are of no interest to us, the lines $\xi_2$ and $\dot{\xi}_3$ are absent from the differential system (7.8). The system (7.8) only contains the tidal perturbations, and the total differential system that we consider for our model is the one derived from the Hamiltonian (4.21), $F_0$, to which we add the tidal perturbations (7.8). We denote it $F: \mathbb{R}^{13} \to \mathbb{R}^{13}$.

## 7.2 Pseudo-equilibria and linearization

### 7.2.1 The eigenvalues

As explained at the beginning of this section, we want to find the equilibria of $F$ and to study the linearized dynamics in their vicinity. However, although $F_0$ has many equilibria (Table 4.1), we show that $F$ has none. The last five lines of Eqs. (7.8), those corresponding
to $\dot{\Gamma}$, $\dot{G}$, $\dot{\omega}_1$, $\dot{\omega}_2$ and $\dot{\omega}_3$, cannot all vanish if $D_j \neq 0$. Indeed, if the planets are all synchronized, that is, if the $\dot{\omega}_j$ all vanish$^3$, then $\dot{\Gamma} < 0$ and

$$ \dot{\Gamma} \approx -\sum_{j=1}^{3} 21 \frac{n_{j,0}}{n_{3,0}} n_{j,0} \frac{m_0}{m_j} \frac{q_j}{Q_j} D_j \propto -\sum_{j=1}^{3} \frac{m_0}{m_j} \frac{q_j}{Q_j} e_j^2. \quad (7.9) $$

However, $F_0$ has no equilibria at $D_j = 0$ (Table 4.1) and does not contribute to these five lines. We thus conclude that $F$ has no equilibria. This result contrasts with the case of co-orbitals without companion, studied in Chap. 6. Indeed, the differential system (6.15), that governs the averaged dynamics of the resonance chain 1 : 1 under tidal dissipation, does have equilibria. This is because the conservative part of the differential system (6.15), derived from the Hamiltonian $\mathcal{H}_K + i\mathcal{H}_P$ in Eq. (6.11), has equilibria at zero eccentricity, that is, at $D_j = 0$, whereas all the equilibria of $F_0$ are at non-zero eccentricity.

We call pseudo-equilibrium of $F$, or pseudo-fixed point of $F$, a point $\mathcal{X} \in \mathbb{R}^{13}$ such that $F(\mathcal{X}) = {}^t(0, 0, 0, 0, 0, 0, 0, \dot{\Gamma}(\mathcal{X}), 0, 0, 0, 0)$. Even though it does not have equilibria, $F$ has pseudo-equilibria, and we find them using an extension of the Newton-Raphson-based algorithm that we developed for Sect. 4.3.2. In the case of the resonance chain $p : p : p+1$, the variable $\Gamma$, that many authors refer to as scaling parameter (Michtchenko et al., 2008; Delisle, 2017; Petit et al., 2020), drifts along the branches of pseudo-equilibria of $F$. In the case of the co-orbital resonance 1 : 1 however, the scaling parameter is$^4 J_2$, and while this quantity drifts at a random point of the phase space, it does not drift along the branches of equilibria of the differential system (6.15).

Equation (7.9) shows that, on a branch of pseudo equilibria of $F$, the parameter $\Gamma$ (and thus the parameter $\delta$, see Eq. (4.35)) drifts at a speed proportionnal to the square of the eccentricities. This means that, with tides, the system travels along the main branch from right to left in Fig. 4.2 (Delisle et al., 2014) much more quickly when $\delta > 0$ than when $\delta < 0$ (due to high eccentricities for positive $\delta$). As the system travels, whether or not it stays close to the main branch or moves away depends on the linear stability of the differential system $F$ in the vicinity of the branch. That is, it depends on the real parts of the eigenvalues of the linear system associated with $F$ (Sect. 2.4.1). Since $\Gamma$ is not constant at the pseudo-fixed points, but drifts at a speed given by Eq. (7.9), computing the eigenvalues of the linearized system makes sense only if $\Gamma$ drifts slowly enough, that is, only if

$$ |\dot{\Gamma}| \ll \max_{1 \le k \le 13} |\mathfrak{Re} \, \lambda_k|, \quad (7.10) $$

where the $\lambda_k$ are the eigenvalues of the linearized system. Indeed, $|\dot{\Gamma}|^{-1}$ is the timescale of evolution of $\Gamma$, while $(\max_{k \le 13} |\mathfrak{Re} \, \lambda_k|)^{-1}$ is the timescale of tidal evolution. When the criterion (7.10) is fulfilled, $\Gamma$ can be considered constant on the timescale of tidal evolution, and the real parts of the linearized system have physical meaning.

Branch 3 always has high eccentricities (see Fig. 4.3) and exists only for $\delta > 5.997$. The drift in $\delta$ towards negative values is fast at high eccentricity (see Eq. (7.9)), and so branch 3 is tidally very unstable and uninteresting to us. Branch 2 has low values of the eccentricities at large $\delta$, but the existence of a hyperbolic zone at $5.55 \le \delta \le 5.80$ (see Table 4.1) makes it uninteresting as well since the drift ensures that this zone is reached. Hence, we limit the study of tidal dissipation to the main branch (branch 1).

<footnote>$^3$Then, $\dot{G}$ also vanishes, since $\dot{G} + \sum_j \alpha_j \dot{\gamma}_j^2 C_j^{-1} \dot{\omega}_j = 0$ by conservation of the total angular momentum.</footnote>

<footnote>$^4$In a resonance chain, the scaling parameter is the momentum associated with the averaged fast angle.</footnote>
<figure><img src="image_22.png" /><figcaption>Fig. 7.1 — Real parts of the eigenvalues of the linearized system associated with $F$ in the vicinity of its main branch of pseudo-equilibria, for $-7 \le \delta \le 1$ (left) and $-6.2 \le \delta \le -5$ (right). The planetary masses are as in Fig. 4.2, the resonance chain is 1 : 1 : 2 and the tidal parameters are those of system 1 in Table 7.1. Only the four eigenvalues (actually eight, counting their complex conjugated) associated with the four degrees of freedom of the conservative system are represented. The five other eigenvalues, associated with the three rotation rates $\omega_j$, with the total angular momentum, and with $\Gamma$, are of no interest. The real parts behave erratically at the 1 : 1 secular resonance between $\nu$ and $\nu_3$ (see Fig. 4.2), in such a way that all of them are negative for $-5.64 \le \delta \le -5.47$, yielding a linear stability in this region (Sect. 2.4.1). In the region $\delta \le 1$ we have $|\dot{\Gamma}|/\Re\lambda \le 0.012$, and the criterion (7.10) is very well respected. Similar figures with other planetary masses are available in Sect. 7.3.4.</figcaption></figure>

### 7.2.2 The linearly stable region

As explained in Sect. 4.3.1, it is not possible to give analytically the equilibria of $F_0$, and so a fortiori, it is not possible to give the equilibria of $F$. Consequently, it is clear that the eigenvalues of the differential system $F$, linearized in the vicinity of the main branch, cannot be given analytically. In particular, the perturbative approach described in Sect. 2.4.2 is of no help here, since we do not have an analytical expression of a diagonal basis of $M_0$ in Eq. (2.82). This is a notable difference with the case of the co-orbital resonance 1 : 1 for which we were able to give analytical expressions of the eigenvalues with tides (Eq. (6.28)).

We work numerically instead, and in Fig. 7.1 we plot the real parts of the eigenvalues of the linearized system associated with $F$, along its main branch of pseudo equilibria, which is a little perturbation of the main branch of equilibria of (4.21). To guarantee that the condition (7.10) is well respected, we limit ourselves to $\delta < 1$. This is not really a restriction since tides ensure that this region is quickly reached. Only the eigenvalues that are the extension of Eqs. (4.38) & (3.32) are plotted. In the absence of a third planet, we showed with Eq. (6.28) that the eigenvalue responsible for the exponential increase in the
<figure><img src="image_23.png" /><figcaption>Fig. 7.2 — Position of the linearly stable region for the resonance chain 1 : 1 : 2 in the plane $(m_3, \delta)$. For every point in this plane the eigenvalues of the linearized system associated with the vector field $F$ are computed, and the point is plotted only if all the real parts are negative (actually, if they all are less than $0.05 \Re \lambda$). The co-orbital masses are $m_1 = m_2 = 10^{-4}$ and the tidal parameters are those of the system 0 in Table 7.1. The position of the linearly stable region weakly depends on $m_1/m_2$ and on the tidal parameters. The colour gives $n_1/n_3$ and shows that the chain stabilizes the dynamics far from the Keplerian resonance (for which $n_1/n_3 = 2$). The dashed yellow line plots the secular 1 : 1 resonance between $\nu$ and $\nu_3$, computed with Eqs. (4.39) and (4.30), respectively. For $m_3 > 18$ ($m_1 + m_2$), the linearly stable region disappears, while for $m_3 < 0.29$ ($m_1 + m_2$), two distinct linearly stable regions exist, whose widths tend to 0 with $m_3$.</figcaption></figure>

libration amplitude of $\xi$, and thus for the destruction of the co-orbital motion, has a real part

$$ \Re \lambda = \frac{9}{2} \eta \frac{m_0}{m_1 + m_2} \left( \frac{m_1 q_2}{m_2 Q_2} + \frac{m_2 q_1}{m_1 Q_1} \right), \qquad (7.11) $$

and so we normalize by this quantity in Fig. 7.1 and Sect. 7.3.4. The region within $-5.64 \le \delta \le -5.47$, at the 1 : 1 secular resonance between $\nu$ and $\nu_3$ (i.e. between the libration frequency of the co-orbitals and the frequency of all the pericentres at the pseudo-equilibria), is such that all the eigenvalues of the system have negative real parts, and we thus expect this region to be linearly stable<fnref n="5" /> (Sect. 2.4.1). We show in Sect. 7.3.2 that this is indeed the case. The linear stability is only temporary though, since the drift in $\delta$ ensures that this region is eventually left. We call it the *linearly stable region* in the remainder of this work. The range in $\delta$ corresponding to the *linearly stable region* strongly depends on the ratio $m_3/(m_1 + m_2)$. In Fig. 7.2 we display its position in the plane $(m_3, \delta)$.

<footnote n="5">⁵Even though it is slightly chaotic far from the main branch, according to the map in Fig. 4.6.</footnote>
## 7.3 Numerical simulations

In this section we investigate the ability of our model to predict the behaviour of a system in the $p : p : p + 1$ resonance under tidal dissipation. We especially check the results drawn in Sect. 7.2.2 on the linearized dynamics in the vicinity of the main branch and the existence of a *linearly stable region*. In Sect. 7.3.4, we investigate the influence of $m_3$ on the co-orbital dynamics.

### 7.3.1 Procedure

We numerically integrate two different sets of equations. The first set, our model, is the differential system $F$, that is, the vector field $F_0$ derived from the conservative Hamiltonian (4.21), to which we add the tidal perturbations given by Eqs. (7.8). The second set is a $n$-body unaveraged and unexpanded model in the regular position-speed coordinates, with the constant-$\Delta t$ model, given by Eq. (C.13) of appendix C.4.

When the third planet is absent, we showed in Sect. 6.2.2 that the relevant parameters to consider to predict the destruction time of a system of two co-orbital planets are (Eq. (6.2.2))

$$ \Omega = \frac{q_1}{Q_1} + \frac{q_2}{Q_2}, \quad x = \frac{m_1}{m_2}, \quad y = \frac{q_2 Q_1}{q_1 Q_2}, \quad \text{and} \quad \iota = \frac{m_1 + m_2}{m_0}, \qquad (7.12) $$

namely the total dissipation rate, the co-orbital mass ratio, the dissipation rate ratio, and the planet-to-star mass ratio. Since we are interested in comparing the lifetimes of the co-orbitals when they are inside the resonance chain 1 : 1 : 2 with their lifetimes when they are alone, we make use of these parameters. As explained in Sect. 7.2.2, we do not have analytical expressions depending on the parameters of the lifetime of the system, and trying to draw a complete picture would require a very large number of simulations. Instead, we are interested in performing a small number of simulations with parameters that we judge interesting. Thus, in order to check the semi-analytical results obtained in Sect. 7.2.2, we only show the evolution of two systems for the chain 1 : 1 : 2, whose parameters are given in Table 7.1. Nevertheless, we performed additional simulations with different choices for the planetary masses and the initial $\delta$. The most interesting ones are presented in Sect. 7.3.4, where we thoroughly discuss the influence of larger or smaller values for $m_3$.

For the systems listed in Table 7.1, the chosen value of $\delta$ is such that the beginning of the simulations is at the rightmost point of the *linearly stable region*; these regions

<table><caption>Table 7.1 — Parameters of the two numerical simulations of Figs. 7.3 & 7.4. In this table, the chosen value of $\delta$ is that of the maximum of the region where all the real parts are negative, ensuring that the system crosses the whole *linearly stable region*. The rotations rates $\omega_j$ are initially synchronized and both sets of parameters verify $m_3/m_0 = 10^{-4}$, $\kappa_2^{(3)} = 0$ and $\bar{a} = 0.02$ AU. System 1 is the system that was used for all the figures in Chap. 4 and Sect. 7.2.</caption><thead><tr><th>#</th><th>Ω</th><th>x</th><th>y</th><th>ι</th><th>δ</th></tr></thead><tbody><tr><td>0</td><td>4 × 10<sup>-12</sup></td><td>1</td><td>1</td><td>2 × 10<sup>-4</sup></td><td>-5.51</td></tr><tr><td>1</td><td>4 × 10<sup>-12</sup></td><td>1/10</td><td>100</td><td>2 × 10<sup>-4</sup></td><td>-5.46</td></tr></tbody></table>
<figure><img src="image_24.png" /></figure>

Fig. 7.3 — Evolution of $\xi$ and the $\sigma_j$ (in arc degrees) as a function of time for system 1 (see Table 7.1) as integrated by the simplified model $F$ (Eqs. (4.21) & (7.8)) (top) and the direct $n$-body simulation (C.13) (bottom). In the direct simulation horseshoe-shaped orbits are reached after 6.42 $\tau_{hs}$, and the co-orbital motion is destroyed shortly after that (chaoticity in Fig. 4.6). The model reaches the horseshoe-shaped orbits at $t = 13.1 \tau_{hs}$. Here, the presence of the chain increases the lifetime of the co-orbitals by a factor of 6.42. The thickness of the lines in the bottom plot (e.g. $\sigma_3$) is due to the short-period oscillations that were averaged out in the model. The grey-shaded area is the *linearly stable region*, that the system leaves at $t = 0.2 \tau_{hs}$.

are $-5.67 \le \delta \le -5.52$ for system 0 and $-5.64 \le \delta \le -5.47$ for system 1 (see Fig. 7.1). These systems are thus expected to be initially very stable until they leave this region (due to the drift in $\delta$, see Sect. 7.2.1).

To integrate the two systems with the simplified model $F$ (Eqs. (4.21) & (7.8)), we find, for the given value of $\delta$ and the planetary masses, the position of the fixed point of $F_0$ along the main branch, and use it as initial condition for the integration, with a shift $\Delta\xi = 0.1^\circ$ in $\xi$, in order to not start exactly at the fixed point. The pseudo-fixed point of $F$ is very close to the fixed point of $F_0$, and we ignore the difference. To integrate the system with the complete $n$-body set of Eqs. (C.13), we find, for the given value of $\delta$ and the planetary masses, the position of the *libration centre* along the main branch with the algorithm described in Sect. 4.2.2, and use it as the initial condition for the integration, again with the shift $\Delta\xi = 0.1^\circ$.

### 7.3.2 Increased co-orbital life inside the resonance chain

When the co-orbital planets are alone, the positivity of $\Re c \lambda$ in Eqs. (7.11) & (6.28) ensures that the system systematically reaches the horseshoe-shaped orbits and is destroyed by close encounters (Fig. 6.3). In this case the time $\tau_{hs}$ needed to reach the horseshoe-shaped orbits is given by Eq. (6.40), while the time $\tau_{dest}$ until destruction is close to $\tau_{hs}$ for a
<figure><img src="image_25.png" /><figcaption>Fig. 7.4 — Evolution of $\xi$ and the $\sigma_j$ (in arc degrees) as a function of time for system 0 as integrated by the simplified model $F$ (Eqs. (4.21) & (7.8)) (top) and by the direct $n$-body simulation (C.13) (bottom). In the direct simulation the horseshoe-shaped orbits are reached after 7.8 $\tau_{hs}$, and the co-orbital motion is destroyed shortly after that (chaoticity in Fig. 4.6). The model reaches horseshoe-shaped orbits at 8.1 $\tau_{hs}$. The grey-shaded area is the *linearly stable region*, that the system leaves at $t = 1.05 \tau_{hs}$.</figcaption></figure>

large range in $\iota$ (Sect. 6.2.3), and we approximate $\tau_{hs} \sim \tau_{dest}$. In order to better compare the co-orbital lifetime inside and outside the chain $p : p : p + 1$, we normalize the time by $\tau_{hs}$ in Figs. 7.3, 7.4, 7.5, and in Sect. 7.3.4.

In Figs. 7.3 and 7.4 we plot the angles $\xi$ and $\sigma_j$ as a function of time. The system spends a large amount of time close to the main branch of equilibria, which allows the co-orbitals to live notably longer with the presence of the third planet. This can be seen from the destruction occurring at a time $t > \tau_{hs}$. When the system crosses the *linearly stable region*, the libration amplitude of $\xi$ decreases instead of increasing exponentially, since the real parts of all the eigenvalues of the linearized system associated with $F$ are negative. When the system leaves this region due to the drift in $\delta$ and the real part of one eigenvalue becomes positive again, the libration amplitude of $\xi$ is much smaller than it was before entering the *linearly stable region*. As a result, the system needs more time to reach large libration amplitudes and settle in horseshoe-shaped orbits, which delays the co-orbital destruction; in other words, crossing the *linearly stable region* while being sufficiently close to the main branch (so that the linear dynamics dominate) guarantees a co-orbital lifetime longer than without the third planet.

In Fig. 7.3, the negative real parts of all the eigenvalues in the *linearly stable region* allow the libration amplitude of $\xi$ to reach values as small as 1.8 arcseconds at $t = 0.57 \tau_{hs}$. This minimum happens after the *linearly stable region* is left since the proper mode associated with the newly positive real part (in yellow in Fig. 7.1) has been completely damped by the *linearly stable region* and some time is needed to pump it noticeably.
<figure><img src="image_26.png" /></figure>

Fig. 7.5 — Value of $e_1$ as a function of time for system 1 in Table 7.1, as integrated by the simplified model *F* (Eqs. (4.21) & (7.8)) (*bottom*), and a schematic representation of the eccentricity vector $e_1e^{i\varpi_1}$, for four particular times (*top*). The schema explains the different stages of the *eccentricity damping stabilization* mechanism. The drift in $\delta$ brings the centre of the circle drawn by $e_1e^{i\varpi_1}$ (equilibrium value of $e_1$ along the main branch) closer to the origin (see Fig. 4.2).

Similarly, in Fig. 7.4, 3.5 arcseconds of libration amplitude are reached at $t = 1.95 \tau_{hs}$. In the bottom panel of Fig. 7.7 from Sect. 7.3.4, we display the important damping of $\xi$. The early augmentation of the libration amplitude of the angles in the bottom plot is due to the fact that in the complete simulation (C.13), the *linearly stable region* is not exactly at the same values of $\delta$ as in the simplified model.

### 7.3.3 The eccentricity damping stabilization mechanism

The *linearly stable region* is not the only reason why co-orbitals in resonant chains can live longer. Another phenomenon, which we refer to as *eccentricity damping stabilization* in the rest of this work, allows the libration amplitude to not cross the separatrix leading to horseshoe-shaped orbits. After an exponential increase in the libration amplitudes of $\xi$ and $\sigma_j$, due to at least one eigenvalue with a strictly positive real part, the amplitudes suddenly decrease and the system returns close to the equilibria. This stabilization of $\xi$, due to eccentricity damping (see Fig. 7.5), can happen several times before horseshoe-shaped orbits are finally reached, and the system is destroyed (see Fig 7.4).

The explanation of the *eccentricity damping stabilization* relies on the behaviour of the eccentricities. In Fig. 7.5, we plot the eccentricity $e_1$ of planet 1 as a function of time, together with a schema of its behaviour in the plane $(e \cos \varpi, e \sin \varpi)$. At time $t = t_1$ the system is still close to the fixed points, hence the quantity $e_1e^{i\varpi_1}$ (or any of the two other eccentricities) describes a circle of small radius. Outside the *linearly stable region*, the eigenvalues of the linearized systems have one positive real part, and as time evolves the radius of the circle grows, while its centre, the equilibrium position of $e_1e^{i\varpi_1}$, gets closer to the origin due to the drift in $\delta$ (the eccentricities are smaller on the main branch at
small $\delta$-values, see Fig. 4.2). At time $t > t_2$ the circle $e_1 e^{i\varpi_1}$ surrounds the origin and keeps growing, as predicted by the eigenvalues, which triggers a jump in the eccentricity. On the one hand, the linearized system predicts that the circle drawn by $e_1 e^{i\varpi_1}$ grows to infinity; on the other hand, tides impose an exponential decay of the eccentricities. Indeed, the first line of Eqs. (7.8) yields (Correia, 2009)

$$ \dot{e}_j = -\frac{e_j}{\tau_j}, \quad \tau_j = \frac{2 m_j Q_j}{21 m_0 q_j n_j^{(0)}}, \qquad (7.13) $$

and so, at time $t = t_3$, the circle reaches its maximum radius, the system is now far from its equilibrium, and tides, through non-linear contributions of the vector field $F$, force the eccentricities to decrease. At $t = t_4$, the system has returned into the vicinity of the main branch, where the libration amplitude of the $\sigma_j$, but also of $\xi$, is small. In Figs. 7.3 and 7.4, several occurrences of the *eccentricity damping stabilization* prevent the angle $\xi$ from reaching the horseshoe-shaped orbits and increase the lifetime of the co-orbital planets.

While the stability induced by the *linearly stable region* comes from linear contributions of the vector field $F$, the *eccentricity damping stabilization* comes from non-linear contributions. This latter mechanism works thanks to a strong coupling between the eccentricities ($D_j; \sigma_j$) and the co-orbital angle ($L; \xi$). In the region $\delta < 0$ (tidally interesting), at most one eigenvalue has a positive real part (see Fig. 7.1), but due to the coupling, it allows an exponential growth of the libration angle $\xi$, as well as the eccentricities, which makes the *eccentricity damping stabilization* possible. When the time $t = t_3$ is reached, the coupling ensures that the eccentricity damping also induces a damping of the libration angle $\xi$, hence the stabilization of the co-orbital motion. In the absence of the third planet, we showed in Sect. 6.2.1 that the eccentricities are uncoupled from the co-orbital angle $\xi$ (the matrix $Q_0 + Q_1$ is block diagonal). This means that the positive real part $\lambda$ (Eqs. (7.11) & (6.28)) associated with ($L; \xi$) does not induce an exponential growth of the eccentricities, which are on the contrary damped to 0 due to the negative real parts of their eigenvalues. In this case, because of the decoupling, even if some other mechanism increases the eccentricities, the eccentricity damping predicted by Eq. (7.13) still occurs, but it does not induce a stabilization of $\xi$.

The occurrence of the *eccentricity damping stabilization* is not systematic. It occurs only if the time $t = t_3$ happens before the co-orbital planets reach horseshoe-shaped orbits. If not, the co-orbital configuration is destroyed before the exponential decrease in the eccentricities can save it. Figures 7.8 & 7.9 from Sect. 7.3.4 give examples of systems where the *eccentricity damping stabilization* fails (systems 3, 4 and 6 in Table 7.2). Deciding whether or not a given system will be saved by the *eccentricity damping stabilization* requires to know the proper modes of the linearized system associated with $F$ and how ($L; \xi$) and the ($D_j; \sigma_j$) are written in the corresponding diagonal basis. Only a numerical work is possible, and we did not undertake it since it is much easier to simply run the corresponding numerical simulation.

If a larger initial $\delta$-value is chosen in these simulations, the system initially has at least one positive real part and moves away from the fixed point at exponential speed. If the *eccentricity damping stabilization* works, or if the initial $\delta$-value is small enough, the *linearly stable region* is reached. However, if the system reaches the *linearly stable region* when it is too far from the equilibria, the non-linear contributions of $F$, combined with the chaotic motion induced by the 1 : 1 secular resonance between $\nu$ and $\nu_3$ (see
<figure><img src="image_27.png" /><figcaption>Fig. 7.6 — Evolution of $\xi$ and the $\sigma_j$ for $m_3 = 8(m_1 + m_2)$ (left) and real parts of the eigenvalues in the vicinity of the main branch (right).</figcaption></figure>

the stability map in Fig. 4.6), can lead to peculiar orbits (e.g. switching between the Lagrangian equilibria $L_4$ and $L_5$, that is, a permutation of the co-orbitals). Entering the *linearly stable region* while still close enough to the equilibria ensures a convergence towards the main branch, and thus an increased stability.

It can be seen in Fig. 4.5 that for $\delta = -5.46$ the system is already far from the exact resonance. For system 1 at $\delta = -5.46$, we have $n_1/n_3 = 2.072$. As time goes by, $\delta$ drifts towards more negative values, and at $t = 6.42 \tau_{hs}$, when it is about to reach horseshoe-shaped orbits and be destroyed, system 1 verifies $\delta = -13.77$ and $n_1/n_3 = 2.186$. Similar considerations are valid for system 0, which means that the system is already outside the resonance, but it is still influenced by the chain. As the system leaves the resonance due to the drift in $\delta$, the coupling between the eccentricities and $\xi$ becomes weaker, meaning that the *eccentricity damping stabilization* ends up failing, although it can still work for values of $n_1/n_3$ as high as 2.6. This prevents the co-orbitals from living forever. Only positive values of $\delta$ allow for $n_1/n_3$ a value close to $(p+1)/p = 2$ (see Fig. 4.5). When a positive value of $\delta$ is chosen at $t=0$, the quick drift in $\delta$ due to the high values of the eccentricities (Eq. (7.9)) forces the system to reach the region $\delta < 0$ on a timescale generally shorter than $1/\Re\lambda$, where it can be saved by the *linearly stable region*. For systems on the main branch, this means that tides favour for the ratio $n_1/n_3$ values above their Keplerian value. This result was shown by Delisle *et al.* (2014) for a two-planet chain and is confirmed by the observations of the Kepler mission, where a large number of exoplanets were discovered with a mean motion ratio slightly higher than $(p+1)/p$ (e.g. Delisle and Laskar, 2014).

This section shows that our averaged model (Eqs. (4.21) & (7.8)) is able to satisfyingly predict the tidal evolution of a resonance chain of the form $p : p : p+1$, at least qualitatively, since a precise quantitative description can only be achieved by running the simulation of the complete set of Eqs. (C.13). This contrasts with the 1 : 1 mean motion resonance
<figure><img src="image_28.png" /><figcaption>Fig. 7.7 — Evolution of $\xi$ and the $\sigma_j$ for $m_3 = 0.29 (m_1 + m_2)$ (left) and real parts of the eigenvalues in the vicinity of the main branch (right).</figcaption></figure>

alone, that we studied in Chap. 6, where the secular model is able to quantitatively predict the outcome of the complete simulations of the unaveraged system with less than 1% relative error (Fig. 6.3 & Table 6.3).

### 7.3.4 Influence of the mass of the third planet

In this section, we present the six most interesting simulations that were not shown earlier in Sect. 7.3. We particularly focus on the influence of the mass $m_3$ on the co-orbital dynamics. All the simulations comply with $m_1 = m_2 = 10^{-4} m_0$, and their tidal parameters are those of system 0 in Table 7.1. We only integrate here the simplified model $F$ (Eqs. (4.21) & (7.8)), defined at the end of Sect. 7.1.2.

For each simulation, the real parts of the eigenvalues of the linearized system associated with $F$ are shown alongside the time evolution of the angles $\xi$ and $\sigma_j$. In the figures of the real parts, a dashed vertical black line shows the starting value of $\delta$ of the corresponding simulation. In the figures of the angles, a grey-shaded area shows the *linearly stable region*, when relevant. Choosing other tidal parameters does not significantly modify the figures shown here, since we normalize the real parts by $\Re\lambda$ (Eqs. (7.11) & (6.28)) and the times by $\tau_{hs}$ (Eq. (6.40)), which is the time to reach horseshoe-shaped orbits (close to the destruction time) in the absence of a third planet (Sect. 6.2.3). The mass of the third planet, $m_3$, and the value of $\delta$ at $t=0$, denoted by $\delta_0$, are the only two varying parameters between the different simulations, and these two quantities are also the axes of Fig. 7.2 that displays the *linearly stable region*. The parameters $m_3$ and $\delta_0$ of the six systems integrated here are given in Table 7.2. For each $\delta_0$ and $m_3$, the simulation starts at the corresponding point of the main branch, with a shift $\Delta\xi = 0.1^\circ$ to $\xi$, in order not to start exactly at the equilibrium.

In Fig. 7.6, we have $m_3 = 8(m_1 + m_2)$, which corresponds to the value yielding the largest *linearly stable region* (Fig. 7.2). However, for this choice of $m_3$, the *linearly stable*
<figure><img src="image_29.png" /><figcaption>Fig. 7.8 — Evolution of $\xi$ and the $\sigma_j$ for $m_3 = (m_1 + m_2)/32$ (left) and real parts of the eigenvalues in the vicinity of the main branch (right).</figcaption></figure>

region is located at larger values of $\delta$, which then drifts quickly (Sect. 7.2.1). Furthermore, Eq. (7.9) shows that the drift in $\delta$ is proportional to $\sum_j m_j^{2/3}$, and so, a rather large $m_3$-value is responsible for a quicker drift. Despite the *linearly stable region* being thickest for this value of $m_3$, system 1 leaves it after less than 0.1 $\tau_{hs}$, much more quickly than in Fig. 7.7, where $m_3$ is smaller. The amplitude of libration of $\xi$, $\Delta\xi = \max(\xi - \xi_{eq})$, reaches 0.07° at its lowest, not significantly smaller than its initial value of 0.1°. However, for this choice of $m_3$, the *eccentricity damping stabilization* is very efficient and occurs seven times, yielding a co-orbital lifetime of 8 $\tau_{hs}$.

In Fig. 7.7, we have $m_3 = 0.29(m_1 + m_2)$, which is the best compromise between the width of the *linearly stable region* and the speed of the drift in $\delta$. System 2 stays in this region for nearly 2 $\tau_{hs}$, 20 times longer than system 1 in Fig. 7.6. The *linearly stable region* is so efficient that at $t = 3.2 \tau_{hs}$, the libration amplitude of $\xi$ reaches $\Delta\xi = 4$ milliarcseconds, gaining a factor $\sim 80000$ from the initial $\Delta\xi = 0.1°$. The system stays close to the main branch for 6 $\tau_{hs}$, which is longer than any other simulation that we performed. We also display in Fig. 7.7 a zoomed-in view of $\xi$ close to the main branch to better visualize the damping of its libration amplitude. At time $t \sim 1 \tau_{hs}$, $\Delta\xi$ temporarily stops decreasing and even slightly increases. Indeed, for this value of $m_3$, the *linearly stable region* splits into two distinct strips (Fig. 7.2), and in the middle of the region, the

<table><caption>Table 7.2 — Parameters of the six systems integrated in Sect. 7.3.4. All six systems verify $m_1 = m_2 = 10^{-4} m_0$, and their tidal parameters are those of system 0 in Table 7.1.</caption><thead><tr><th></th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th></tr></thead><tbody><tr><td>$m_3 / (m_1 + m_2)$</td><td>8</td><td>0.29</td><td>1/32</td><td>1/32</td><td>22</td><td>22</td></tr><tr><td>$\delta_0$</td><td>-1.466</td><td>-6.66</td><td>-10.735</td><td>-10.6</td><td>-1.3</td><td>-7.4</td></tr></tbody></table>
<figure><img src="image_30.png" /><figcaption>Fig. 7.9 — Evolution of $\xi$ and the $\sigma_j$ for $m_3 = 22(m_1 + m_2)$ (left) and real parts of the eigenvalues in the vicinity of the main branch (right).</figcaption></figure>

real part of $\lambda_2$ (the blue eigenvalue in Fig. 7.7) temporarily becomes slightly positive and reaches 0.067 $\Re\lambda$. This is in perfect agreement with the slight increase observed for $\Delta\xi$. The *eccentricity damping stabilization* is not as efficient as it was for system 1, but still allows the system to reach 8.5 $\tau_{hs}$ before the destruction of the co-orbital motion.

In Fig. 7.8, we have $m_3 = (m_1 + m_2)/32$. With such a low mass for the third planet, the two *linearly stable regions* are extremely narrow, and the linear stability occurs only if the system is exactly at the 1 : 1 secular resonance between $\nu$, the libration frequency of $\xi$ (Eq. (3.32)), and $\nu_3$, the precession frequency of the $\varpi_j$ (Eq. (4.30)). With such a tiny difference on $\delta_0$, the initial conditions between system 3 and system 4 are very close, but the outcomes are completely different, as system 3 lives three times longer than system 4. Outside the 1 : 1 secular resonance between $\nu$ and $\nu_3$, the unique positive real part is equal to $\Re\lambda$ and the horseshoe-shaped orbits are reached at $\tau_{hs}$ exactly, as in the absence of the third planet. This is consistent with the small value of $m_3$. The *eccentricity damping stabilization* fails for both systems.

In Fig. 7.9, the third planet has a mass $m_3 = 22(m_1 + m_2)$ and the *linearly stable region* does not exist (Fig. 7.2). We perform a simulation at the 1 : 1 secular resonance between $\nu$ and $\nu_3$, where it would have been if it existed. In this region the positive real parts are not greater than $\Re\lambda$, and with the action of the *eccentricity damping stabilization*, system 5 is destroyed after more than 3 $\tau_{hs}$. However, system 6 is very far from the Keplerian resonance and $\lambda_4$ (the yellow eigenvalue in Fig. 7.9) has a real part $4\Re\lambda$. This means that this system reaches the horseshoe-shaped orbits after a time $\tau_{hs}/4$ only. Furthermore, since $(L, \xi)$ and $(D_j, \sigma_j)$ are weakly coupled far from the Keplerian resonance, the *eccentricity damping stabilization* fails, leading to a very premature destruction of the co-orbital motion (the co-orbitals would have lived four times longer without the third planet). Choosing a smaller $\delta$ or a larger $m_3$ leads to an even quicker destruction.

In this subsection, we showed that the *linearly stable region* is most efficient in
preventing the destruction of the co-orbitals at low to intermediate values of $m_3$, while the *eccentricity damping stabilization* works best at moderately large values of $m_3$, especially if the system is not too far from the Keplerian resonance ($\delta$ not too negative). Generally, a first-order mean-motion resonance with an outward third planet either increases or does not change the co-orbital lifetime. The third planet can nevertheless lead to a premature destruction if it is much more massive than the co-orbitals and far from the mean-motion resonance.

## 7.4 Discussion and conclusion

In this chapter we pushed further the study of the resonance chain $p : p : p + 1$, started in Chap. 4, where an outer third planet is in a first-order mean motion resonance with the pair of co-orbital planets. We took into account dissipation due to tides and we built the pseudo-Hamiltonian of the problem, from which we derived the complete expression of the averaged equations of motion, at second order in eccentricity.

In the conservative case, the secular $1:1$ resonance between the libration frequency of the co-orbital angle and the precession frequency of the pericentres tends to destabilize the chain by chaoticity. When tidal dissipation is involved though, we showed that this secular resonance greatly stabilizes the chain by making negative the real part of the eigenvalues in the vicinity of the main branch of equilibria. This negativity of the real part is responsible for the existence of a value for $n_1/n_3$, close to the secular resonance, around which the libration amplitude of the co-orbital angle $\lambda_1 - \lambda_2$ is damped to zero. As a consequence, the lifetime of the co-orbital pair is noticeably longer with the presence of the third planet. For the chain $1:1:2$, this linear stability generally occurs near $n_1/n_3 \sim 2.07$. Another stabilization mechanism, due to eccentricity damping, is presented and explained. The model that we built is able to predict the position of the *libration centres* of the complete system. When tides are involved, the model reliably gives the qualitative behaviour of the system and, to a certain extent, its quantitative behaviour. This work shows that when tidal dissipation is included, co-orbital systems are more stable if they are inside a resonance chain of the form $p : p : p + 1$, rather than alone, which increases the chances of a future detection of a co-orbital pair of exoplanets.

One important contribution of the work done in Chaps. 4 & 7 is the discovery of the $1:1$ secular resonance between the libration of $\lambda_1 - \lambda_2$ and the precession of the pericentres, as well as the inherent dynamical consequences (Figs. 4.6 & 7.1). *Libration centres* correspond to periodic orbits of the unaveraged problem that generalize the equilibria of the averaged model. They are such that all the pericentres precess at the same frequency in their vicinity, a result that holds true for every resonance chain of any number of planets. We thus expect, for other resonance chains, the existence of similar secular resonances, for example between $\xi = \lambda_1 - 4\lambda_2 + 3\lambda_3$ and the $\varpi_j$ for the chain $1:2:3$.

We proved that tides are responsible for the system traveling towards the left of the main branch in Fig. 4.2. In our model, we only considered the secular contributions due to the resonance $p : p + 1$ (Eqs. (4.19) & (4.20)), which means that the main branch is asymptotic to the horizontal axis and goes up to $\delta = -\infty$. However, Malhotra (2022) showed that the main branch of the resonance $p : p + 1$ is connected to the branch 2
(the other elliptic branch of Fig. 4.2) of the resonance $p-1 : p$. Under the action of the
*linearly stable region* (Sect. 7.2.2) and the *eccentricity damping stabilization* (Sect. 7.3.3),
it is possible, with some choices of the parameters, that before the destruction of the
*co-orbital configuration*, tides allow the system to reach the branch 2 of the $p : p : p+1$
resonance, while traveling from the main branch of the $p+1 : p+1 : p+2$ resonance. We
*did not investigate this scenario here*.
Chapter 8

Conclusions

8.1 Overview

In this thesis manuscript, we studied the dynamics of the co-orbital motion with a
particular focus on the secular perturbations due to tides and mean-motion resonance
with a third planet. As our work makes heavy use of the mathematical theory of celestial
mechanics, we recalled the main parts of this theory in Chap. 2. We rederived in Chap. 3
the important analytical results describing the point-mass 1 : 1 mean-motion resonance.
In particular, we studied there the circular and eccentric dynamics of the co-orbital pair,
which allowed us to define and explain notions that were essential later in the manuscript.
Chapter 3 also contains work original to this thesis and in particular, we showed how
perturbations due to general relativity cannot perturb the co-orbital motion in the same
way as tides do.

We considered once again the point-mass approximation in Chap. 4, where we
switched from a three-body problem to a four-body problem by adding a third planet
to the planetary system. We chose to put the third planet in a first order mean-motion
resonance with the co-orbital pair in order to study the dynamics of a co-orbital pair
perturbed by a resonance chain. As the third planet was exterior to the co-orbital pair,
we studied a resonance chain of the form *p*: *p*: *p* + 1, for any integer *p*. Chapter 4 showed
us that on certain aspects, the dynamics of the point-mass 1 : 1 resonance is similarly
to that of the *p*: *p*: *p* + 1 resonance chain, as simplifications allow the dynamics of the
librating angle $\lambda_1 - \lambda_2$ to be uncoupled from the rest of the dynamics. Nevertheless,
we also discovered specificities inherent in the chain *p*: *p*: *p* + 1, like the existence of
a secular 1 : 1 resonance between the precession frequency of the pericentres and the
libration frequency of the angle $\lambda_1 - \lambda_2$. The vicinity of this secular resonance can make
the trajectories of the system chaotic and we suspect that it exists for other resonance
chains.

In order to properly deal with tidal dissipation, we dedicated the entire Chap. 5 to
revise the main theory on tides. We redefined there important tidal quantities like the
Love numbers and the quality factor, and we studied the tidal models associated with
different rheologies. In particular, we showed how the choice of the constant-$\Delta t$ tidal model allows the equations of motions to be derived from a pseudo-Hamiltonian, as if tides were a conservative perturbation.

In Chaps. 6 & 7, we included tides in the analytical models built in Chaps. 3 & 4, respectively. In Chap. 6, we linearized the equations of motion in the vicinity of the Lagrangian stable equilibria and an analysis of the real parts of the eigenvalues showed us that they are made unstable by tides, whatever are the masses, radii and tidal parameters of the co-orbitals. The instability is due to an exponential increase of the libration amplitude of the angle $\lambda_1 - \lambda_2$, that eventually leads to close encounters between the co-orbitals and the destruction of the system. The amplitude of the real parts gave us informations on the timescale of instability and we were able to give an analytical expression for the lifetime of the co-orbital pair in Eq. (6.40). We explain in Sect. 6.2.3 how this work can be applied to an intelligent search for co-orbital pairs of exoplanets. The analysis of the real part of the eigenvalues is not limited to the obtention of the timescale of destruction but also shows how the eccentric proper modes Lagrange and anti-Lagrange, discovered numerically by Giuppone et al. (2010) and analytically by Robutel and Pousse (2013) are perturbed by tides. We reached the conclusion that the most stable co-orbital pairs settle in the anti-Lagrange configuration. We showed how for a system star–planet–planet, it is enough to consider tides raised on the co-orbitals only, but we also considered in Sect. 6.4 tides raised on the central body to treat the case of a system planet–satellite–satellite made up of small satellites.

In Chap. 7, we included tides in the $p:p:p+1$ resonance chain. We linearized the equations of motions in the vicinity of the equilibria of the system, and we showed that close to the secular 1:1 resonance between the precession of the pericentres and the libration of $\lambda_1 - \lambda_2$, all the real part are negative and the equilibria are made asymptotically stable by tides. For the chain 1:1:2, such secular resonance occurs near $n_3/n_1 \approx 2.07$, outside the Keplerian resonance, and in this region, the libration amplitude of $\lambda_1 - \lambda_2$, instead of increasing exponentially, decreases exponentially. For this reason, the third planet has the ability to greatly extend the lifetime of the co-orbital pair, although it cannot save it from the destruction, as the system always leaves this region of asymptotic stability. We thoroughly studied the influence of the mass of the third planet and we showed that it generally increases the lifetime of the co-orbital pair, although in some few cases, the third planet can lead to a premature destruction of the co-orbital pair.

## 8.2 Consequences

This work, although it shows the systematic destruction of the co-orbital motion, does not entirely explain the absence of detection of co-orbital pair of exoplanets. Indeed, the timescale of destruction is heavily dependent on the semimajor axis (Eq. (6.36)), and co-orbital exoplanets located more than $\sim 0.08$ AU away from their host star can be considered immune to tides. As of September 2, 2022, 52.4% of the 5159 detected exoplanets had a semimajor axis larger than that<fnref n="1" /> and were undergoing from their host star negligible tidal effects. It is thus likely that the challenge that a co-orbital detection represents constitutes the main explanation to the absence of co-orbital exoplanet discovery.

<footnote n="1">¹exoplanet.eu, considering a sample of 2841 exoplanets with known semimajor axis.</footnote>
The disruption due to tidal dissipation constitutes another satisfying explanation. It is also possible that the current formation models overestimate their formation, and co-orbital exoplanets could be rarer than we believe. The results of this work hence cannot be interpreted as a reason to abandon the search for co-orbital exoplanet, and further endeavours will most likely lead to their detection.

In order to prove the unavoidable destruction of the co-orbital motion, this work assumes a very simple tidal model. However, any tidal model leads to the loss of orbital energy, and according to the results of Moeckel (2017) (last paragraph of Sect. 6.2.1), considering a more subtle model could change the timescale of destruction, but it will not save the co-orbital pair. Three more assumptions of our model, namely the expansion in semi-major axis (Eqs. (3.9) & (3.11)), the expansion in eccentricities (Eq. (3.29)), and the averaging process (Eq. (3.17)) were shown to be perfectly valid in Sect. 6.3. The fifth and last approximation that we made is the planar state of the system. In other words, we assumed that a small mutual inclination between the co-orbitals is not pumped by tides like the libration amplitude of the angle $λ_1 – λ_2$. While we did not check this hypothesis, Dobrovolskis and Lissauer (2022) showed that, at least in the restricted case, the mutual inclination between the co-orbital bodies is damped by tides, which validates our last assumption. In the absence of more planets in the planetary system, it is therefore legitimate to consider that our work closes the question of tides raised on co-orbital planets.

## 8.3 Perspectives

While the question of tides in the 1 : 1 mean motion resonance is closed, the work performed in this thesis only briefly explores the case of larger resonance chains. We restricted ourselves to the chain $p : p : p + 1$, for which we showed that whereas the co-orbital lifetime was generally enhanced by the third planet, the destruction was still unavoidable. It is possible that other resonance chains could make the co-orbital motion stable over arbitrarily large timescales, or even indefinitely. In particular, we did not investigate the dynamics of the chain $p : p + 1 : p + 1$, for which the third planet is innermost. It is not absurd to believe that it could have dynamical features completely different from those of the chain $p : p : p + 1$, since such phenomenon has already been observed, even for a two-planet system. Indeed, Beaugé *et al.* (2006) showed that in the 1 : 2 mean motion resonance, the second fundamental model of resonance of Henrard and Lemaitre (1983) ceases to be valid at large eccentricities if the most massive planet is innermost, while it stays valid if the most massive planet is outermost. Something similar is not to be excluded between the chains $p : p : p + 1$ and $p : p + 1 : p + 1$, which is a perspective worth exploring.

Another interesting direction to explore is the scenario mentioned in the last paragraph of Sect. 7.4. As proven by Malhotra (2022), there exist continuous connections between the branches of equilibria of consecutive first-order mean-motion resonances. It would be interesting to build a secular model that takes into account effects of both the $p : p : p + 1$ and $p + 1 : p + 1 : p + 2$ resonance chains, to check if tides, through *linearly stable regions* (Sect. 7.2.2) and *eccentricity damping stabilization mechanisms* (Sect. 7.3.3), could be responsible for the system traveling from the resonance $p + 1 : p + 1 : p + 2$ to the resonance
$p : p : p + 1.$

The case of a system planet-satellite-satellite requires tides raised on the central planet to be taken into account (Sect. 6.4.1). We studied these types of tides in Sect. 6.4, but we limited ourselves to the obtention of analytical expressions of the timescales of evolution. We did not check these analytical results with numerical simulations. Furthermore, for such a system, other effects, like the oblateness of the central stars, or the presence of a non-resonant, much more massive satellite (e.g. the system Saturn-Janus-Epimetheus-Titan) need to be taken into account in order to build a reliable model. We plan to look in this direction in the future.

Finally, the subject of resonance chains is very broad, and for four or more planets, their dynamics are poorly understood, even in the point-mass case. Indeed, adding planets to the chain quickly increases the number of degrees of freedom of the system, and analytical models soon become untractable. Choosing a resonance chain that features a co-orbital pair is a way to reduce the complexity of the model by means of uncoupling simplifications like those explained in Sect. 4.3.1. Furthermore, the presence of a co-orbital pair induces the existence of a fundamental frequency independent of the distance to the Keplerian resonance, namely the libration frequency of the co-orbital angle. The precession frequency of the pericentres, however, depends on the distance to the Keplerian resonance. The presence of a co-orbital pair hence increases the chance of encountering a 1 : 1 secular resonance at a certain distance from the Keplerian resonance, like the one that exists at $n_1/n_3 \approx 2.07$ for the chain 1 : 1 : 2. As an example, the four-planet chain $p : p : p+1 : p+1$ could be comparatively easier to study than other four-planet resonance chains, while being very rich due to two 1 : 1 secular resonances, the first one between $\lambda_1 - \lambda_2$ and the pericentres, and the second one between $\lambda_3 - \lambda_4$ and the pericentres.
# Appendix A

## Spherical harmonics and demonstration of Eq. (5.8)

### A.1 Framework

In this appendix, in order to have more compact notations, we denote

$$ \int_{\mathbb{S}} := \int_{0}^{\pi} \int_{0}^{2\pi} \sin \theta d\varphi d\theta, \quad \text{or similarly} \quad \int_{\mathbb{S}'} := \int_{0}^{\pi} \int_{0}^{2\pi} \sin \theta' d\varphi' d\theta'. \qquad (A.1) $$

We use the spherical harmonics defined by

$$ Y_{lm}(\theta, \varphi) = P_{lm}(\cos \theta) e^{im\varphi}, \qquad (A.2) $$

where the $P_{lm}$ are the associated Legendre polynomials, related to the Legendre polynomials by

$$ P_{lm}(X) = (-1)^m (1 - X^2)^{m/2} \frac{d^m}{dX^m} P_l(X). \qquad (A.3) $$

The Legendre polynomials can be defined in many ways, for example as the coefficients of the serie expansion (5.1), but their most compact definition is Rodrigues' formula

$$ P_l(X) = \frac{1}{2^l l!} \frac{d^l}{dX^l} (X^2 - 1)^l. \qquad (A.4) $$

For two functions $f(\theta, \varphi)$ and $g(\theta, \varphi)$ on the sphere $\mathbb{S}$, we define the dot product

$$ f \cdot g = \int_{\mathbb{S}} f(\theta, \varphi) \bar{g}(\theta, \varphi), \qquad (A.5) $$

where the upper bar is the complex conjugated. With this dot product, the spherical harmonics are orthogonal, that is

$$ Y_{lm} \cdot Y_{l'm'} = c_{lm}^{-1} \delta_{ll'} \delta_{mm'}, \quad \text{with} \quad c_{lm} = \frac{(2l+1)(l-m)!}{4\pi(l+m)!}. \qquad (A.6) $$
The spherical harmonics defined above form an orthogonal basis of the Hilbert space $L_\mathcal{C}^2(\mathbb{S})$ of square integrable functions from $\mathbb{S}$ to $\mathbb{C}$. A given function $f: \mathbb{S} \to \mathbb{C}$ can be decomposed in a unique<fnref n="1" /> way as a linear combination of the $Y_{lm}$

$$f(\theta, \varphi) = \sum_{l=0}^{+\infty} \sum_{m=-l}^{l} f_{lm} Y_{lm}(\theta, \varphi), \qquad (A.7)$$

where the coefficients $f_{lm}$ are given by the projection<fnref n="2" /> of $f$ onto $Y_{lm}$

$$f_{lm} = c_{lm} f \cdot Y_{lm} = c_{lm} \int_{\mathbb{S}} f(\theta, \varphi) \bar{Y}_{lm}(\theta, \varphi). \qquad (A.8)$$

## A.2 Laplace equation in spherical coordinates

The Laplace equation reads $\Delta = 0$. Using the spherical coordinates $\mathbf{r} = (r, \theta, \varphi)$, where $\theta$ and $\varphi$ are the colatitude and longitude, the Laplace operator reads

$$\Delta = \frac{1}{r^2} \frac{\partial}{\partial r} \left( r^2 \frac{\partial}{\partial r} \right) + \frac{1}{r^2 \sin \theta} \frac{\partial}{\partial \theta} \left( \sin \theta \frac{\partial}{\partial \theta} \right) + \frac{1}{r^2 \sin^2 \theta} \frac{\partial^2}{\partial \varphi^2}. \qquad (A.9)$$

Assuming the separation of variables, that is

$$F(r, \theta, \varphi) = f(r)g(\theta)h(\varphi), \qquad (A.10)$$

the equation $\Delta F = 0$ yields

$$f(r) = r^{\beta}, \quad g(\theta) = P_{lm}(\cos \theta), \quad h(\varphi) = e^{\pm im\varphi}, \qquad (A.11)$$

where $l$ is a positive integer, $\beta$ is either $l$ or $-l-1$ and $-l \le m \le l$.

## A.3 Demonstration of Eq. (5.8)

Here we show Eq. (5.8) from Eqs. (5.6) & (5.7). As in Chap. 5, $\sigma$ stands for $(\theta, \varphi)$ and $\mathbb{S}$ denotes the sphere of radius $R$ centered on $\mathcal{O}$, the barycentre of the body $\mathbb{B}$. The Love number $k(R, \sigma, \sigma', t-t')$, appearing in Eq. (5.6), is decomposed over the spherical harmonics

$$\begin{aligned} k(R, \sigma, \sigma', t-t') &= \sum_{l=0}^{+\infty} \sum_{m=-l}^{l} \mathfrak{t}_l^m(R, \sigma', t-t') Y_{lm}(\sigma) \\ &= \sum_{l=0}^{+\infty} \sum_{m=-l}^{l} \left( \sum_{l'=0}^{+\infty} \sum_{m'=-l'}^{l'} \mathfrak{t}_{ll'}^{mm'}(R, t-t') \bar{Y}_{l'm'}(\sigma') \right) Y_{lm}(\sigma). \end{aligned} \qquad (A.12)$$

<footnote n="1">¹$\bar{Y}_{lm}$ can be written in Eq. (A.7) instead of $Y_{lm}$, in which case $f_{lm} = c_{lm} f \cdot \bar{Y}_{lm}$.</footnote>

<footnote n="2">²The coefficient $c_{lm}$ is due to the fact that the $Y_{lm}$ do not form an orthonormal basis of $L_\mathcal{C}^2(\mathbb{S})$, but rather an orthogonal basis.</footnote>
This allows Eq. (5.6) to be rewritten as

$$ 
\begin{aligned}
V(R, \sigma, t) &= \int_{-\infty}^{t} \int_{S'} \sum_{l=0}^{+\infty} \sum_{m=-l}^{l} \left( \sum_{l'=0}^{+\infty} \sum_{m'=-l'}^{l'} \mathfrak{e}_{ll'}^{mm'}(R, t-t') \bar{Y}_{l'm'}(\sigma') \right) Y_{lm}(\sigma) W(R, \sigma', t') dt' \\
&= \sum_{l=0}^{+\infty} \sum_{m=-l}^{l} \int_{-\infty}^{t} \int_{S'} \sum_{l'=0}^{+\infty} \sum_{m'=-l'}^{l'} \mathfrak{e}_{ll'}^{mm'}(R, t-t') \bar{Y}_{l'm'}(\sigma') W(R, \sigma', t') dt' Y_{lm}(\sigma).
\end{aligned} 
\quad (A.13) $$

A term-wise identification between Eq. (5.7) and Eq. (A.13) yields, by uniqueness of the decomposition

$$ V_{lm}(R, t) = \int_{-\infty}^{t} \int_{S'} \sum_{l'=0}^{+\infty} \sum_{m'=-l'}^{l'} \mathfrak{e}_{ll'}^{mm'}(R, t-t') \bar{Y}_{l'm'}(\sigma') W(R, \sigma', t') dt'. \quad (A.14) $$

If $\mathbb{B}$ is isotropic at rest, the coefficients $\mathfrak{e}_{ll'}^{mm'}$ do not depend on the reference frame attached to $\mathbb{B}$, and there exist functions $\mathfrak{e}_l(R, t)$ such that (see Boué et al., 2019)

$$ \mathfrak{e}_{ll'}^{mm'}(R, t-t') = \delta_{ll'} \delta_{mm'} \mathfrak{e}_l(R, t-t'). \quad (A.15) $$

Equation (A.14) hence reduces to

$$ V_{lm}(R, t) = \int_{-\infty}^{t} \mathfrak{e}_l(R, t-t') \int_{S'} \bar{Y}_{lm}(\sigma') W(R, \sigma', t') dt'. \quad (A.16) $$

Substituting Eq. (A.8) into the last equation and defining $k_{lm}(R, t) = c_{lm}^{-1} \mathfrak{e}_l(R, t)$ yields

$$ V_{lm}(R, t) = \int_{-\infty}^{t} k_{lm}(R, t-t') W_{lm}(R, t') dt', \quad (A.17) $$

which is Eq. (5.8).
# Appendix B

## Coefficients of $\mathcal{H}^{(4)}$, $U_t^{(j)}$, $\mathcal{H}_{j,3}^{(1)}$ and $\mathcal{H}_{j,3}^{(2)}$

### B.1 Coefficients of $\mathcal{H}^{(4)}$

The coefficients appearing in the Hamiltonian $\mathcal{H}^{(4)}$ defined in Eq. (3.29) are given by<fnref n="1" />

$$ 
\begin{aligned}
D_h &= \frac{7}{16} \cos \xi + \frac{1}{4\Delta^9} \left( -\frac{3951}{32} + 115 \cos \xi + \frac{293}{8} \cos 2\xi - 27 \cos 3\xi - \frac{37}{32} \cos 4\xi \right), \\
G_h &= \cos \xi + \Delta^{-9} \left( -\frac{4491}{32} + 139 \cos \xi + \frac{233}{8} \cos 2\xi - 27 \cos 3\xi - \frac{25}{32} \cos 4\xi \right), \\
E_h &= \frac{1}{32} \left( e^{-i\xi} + 81e^{-3i\xi} \right) + \frac{e^{-6i\xi}}{32\Delta^9} P_E(e^{i\xi}), \\
F_h &= -\frac{7}{4} e^{2i\xi} + \frac{e^{-3i\xi}}{4\Delta^9} P_F(e^{i\xi}),
\end{aligned}
\quad (B.1) $$

where $P_E$ and $P_F$ are polynomials given by

$$ 
\begin{aligned}
P_E(X) &= -\frac{9}{8} + 15X - \frac{349}{2}X^2 + 171X^3 + \frac{2889}{4}X^4 - 1571X^5 + \frac{2007}{2}X^6 - 87X^7 - \frac{625}{8}X^8, \\
P_F(X) &= \frac{207}{32} + \frac{303}{8}X - \frac{577}{4}X^2 + \frac{603}{8}X^3 + \frac{2511}{16}X^4 - \frac{1475}{8}X^5 + 45X^6 + \frac{57}{8}X^7 - \frac{5}{32}X^8.
\end{aligned}
$$

### B.2 Coefficients of $\mathcal{H}_{j,3}^{(1)}$ and $\mathcal{H}_{j,3}^{(2)}$

We give here the expressions of the coefficients appearing in Eqs. (4.19) & (4.20) of Sect. 4.1.3. They depend on the Laplace coefficients $b_{n/2}^m(\alpha)$ (defined by Eq. (2.75)), and to improve readability, we denote $b_n^m = b_{n/2}^m(\alpha)$, where $\alpha = \bar{a}/a_{3,0}$. For the resonance 1 : 1 : 2

<footnote n="1">¹Recall that $\Delta = \sqrt{2 - 2 \cos \xi}$.</footnote>
we have<fnref n="2" />

$$ 
\begin{aligned}
C_{1,1}^{(1)} &= -\alpha b_3^1 \left( \frac{7}{6} + \frac{2}{3}\alpha^{-2} + \frac{5}{3}\alpha^2 \right) + b_3^0 \left( 1 + \frac{5}{2}\alpha^2 \right) \approx 1.1904937, \\
C_{1,2}^{(1)} &= b_3^1 \left( 1 + \frac{3}{2}\alpha^2 \right) - \frac{5}{2}\alpha b_3^0 + \frac{1}{\sqrt{\alpha}} \approx -0.4283898,
\end{aligned}
\quad (B.2) $$

for the first degree in eccentricity and

$$ 
\begin{aligned}
C_{1,1}^{(2)} &= \alpha b_3^1 \left( \frac{263}{168} + \frac{16}{35}\alpha^{-4} + \frac{89}{105}\alpha^{-2} + \frac{341}{105}\alpha^2 + \frac{184}{35}\alpha^4 \right) \\
&\quad - b_3^0 \left( \frac{71}{70} + \frac{24}{35}\alpha^{-2} + \frac{67}{35}\alpha^2 + \frac{276}{35}\alpha^4 \right) \approx -1.6957266, \\
C_{1,2}^{(2)} &= \alpha b_3^1 \left( \frac{65}{24} + \frac{4}{3}\alpha^{-2} + \frac{13}{3}\alpha^2 \right) - b_3^0 \left( 2 + \frac{13}{2}\alpha^2 \right) \approx -3.5937942, \\
C_{1,3}^{(2)} &= -b_3^1 \left( \frac{29}{10} + \frac{8}{5}\alpha^{-2} + \frac{59}{10}\alpha^2 + \frac{48}{5}\alpha^4 \right) + \alpha b_3^0 \left( \frac{69}{20} + \frac{12}{5}\alpha^{-2} + \frac{72}{5}\alpha^2 \right) \approx 4.9668470, \\
C_{1,4}^{(2)} &= -\frac{1}{8}\alpha b_3^1 \approx -0.3876274, \\
C_{1,5}^{(2)} &= \frac{1}{2}b_3^1 \left( 1 + \alpha^2 \right) - \frac{3}{4}\alpha b_3^0 \approx 0.5756950,
\end{aligned}
\quad (B.3) $$

for the second degree.

## B.3 Coefficients of $\mathcal{U}_t^{(j)}$

We give here the complete expression of the fourth order in eccentricity of the tidal perturbation $\mathcal{U}_t^{(j)}$ to the Hamiltonian, introduced by Eq. (6.8). The term $\mathfrak{D}_4^{(j)}$ appearing in this equation reads

$$ 
\begin{aligned}
\mathfrak{D}_4^{(j)} = D_t^{(j)} \left( \mathcal{R}_j^{-2} X_j^2 \bar{X}_j^2 + \mathcal{R}_j^{*-2} X_j^{*2} \bar{X}_j^{*2} \right) + \left( \mathcal{R}_j \mathcal{R}_j^* \right)^{-1} \left( E_t^{(j)} X_j^2 \bar{X}_j^{*2} + \bar{E}_t^{(j)} X_j^{*2} \bar{X}_j^2 \right) \\
+ \left( \mathcal{R}_j \mathcal{R}_j^* \right)^{-1} G_t^{(j)} X_j X_j^* \bar{X}_j \bar{X}_j^* + \mathcal{R}_j^{-3/2} \mathcal{R}_j^{*-1/2} \left( F_t^{(j)} X_j X_j^* \bar{X}_j^2 + \bar{F}_t^{(j)} X_j^2 \bar{X}_j \bar{X}_j^* \right) \\
+ \mathcal{R}_j^{-1/2} \mathcal{R}_j^{*-3/2} \left( F_t^{(j)} X_j^{*2} \bar{X}_j \bar{X}_j^* + \bar{F}_t^{(j)} X_j X_j^* \bar{X}_j^{*2} \right),
\end{aligned}
\quad (B.4) $$

where

$$ 
\begin{aligned}
D_t^{(j)} &= \frac{3}{8} + \frac{69}{64} \cos [2 (\Delta\lambda_j - \Delta\theta_j)], \\
G_t^{(j)} &= \frac{9}{16} + \frac{75}{16} \cos [2 (\Delta\lambda_j - \Delta\theta_j)], \\
E_t^{(j)} &= \frac{81}{64} e^{-2i\Delta\lambda_j} + \frac{867}{32} e^{-2i(2\Delta\lambda_j - \Delta\theta_j)}, \\
F_t^{(j)} &= \frac{9}{16} e^{i\Delta\lambda_j} - \frac{3}{128} e^{-i(\Delta\lambda_j - 2\Delta\theta_j)} - \frac{1365}{128} e^{i(3\Delta\lambda_j - 2\Delta\theta_j)},
\end{aligned}
\quad (B.5) $$

with $\Delta\lambda_j = \lambda_j - \lambda_j^*$ and $\Delta\theta_j = \theta_j - \theta_j^*$.

<footnote n="2">²Note that these analytical expressions are only valid for $p = 1$, that is, for the resonance chain 1 : 1 : 2.</footnote>
Matrices and differential systems

C.1 Simplified differential system

We give here the expression of the differential system used to obtain the analytical results of Sect. 4.3.1. It is derived from Eq. (2.18) with the Hamiltonian $\mathcal{H}_K + \mathcal{H}^{(0)} + \mathcal{H}^{(1)}$, that is, the Hamiltonian (4.21) truncated at first order in eccentricity, after the degree of freedom ($L$; $\xi$) have been uncoupled from the degrees of freedom ($D_j$; $\sigma_j$), as explained in Sect. 4.3.1. We have

$$
\begin{align*}
\dot{D}_1 &= -\frac{C_{p,1}^{(1)}\sqrt{2C_1D_1}m_1n_{3,0}}{m_0C_3} \sin\left(p\frac{\pi}{3} - \sigma_1\right), \\
\dot{D}_2 &= -\frac{C_{p,1}^{(1)}\sqrt{2C_2D_2}m_2n_{3,0}}{m_0C_3} \sin\sigma_2, \\
\dot{D}_3 &= -\frac{\sqrt{2C_3D_3}C_{p,2}^{(1)}n_{3,0}}{m_0C_3} \left(m_1\sin\left(p\frac{\pi}{3} - \sigma_3\right) - m_2\sin\sigma_3\right), \\
\dot{\Delta}L &= \eta \frac{m_2\sin\xi}{m_0C_1} \left(1 - \frac{1}{\Delta^3}\right), \tag{C.1} \\
\dot{\sigma}_1 &= \frac{\partial \mathcal{H}_K}{\partial \Delta\Upsilon} + \frac{C_1C_{p,1}^{(1)}m_1n_{3,0}}{C_3\sqrt{2C_1D_1}m_0} \cos\left(p\frac{\pi}{3} - \sigma_1\right), \\
\dot{\sigma}_2 &= \frac{\partial \mathcal{H}_K}{\partial \Delta\Upsilon} + \frac{C_2C_{p,1}^{(1)}m_2n_{3,0}}{C_3\sqrt{2C_2D_2}m_0} \cos\sigma_2, \\
\dot{\sigma}_3 &= \frac{\partial \mathcal{H}_K}{\partial \Delta\Upsilon} + \frac{C_{p,2}^{(1)}n_{3,0}}{m_0\sqrt{2C_3D_3}} \left(m_1\cos\left(p\frac{\pi}{3} - \sigma_3\right) + m_2\cos\sigma_3\right), \\
\dot{\xi} &= -3\eta(C_1 + C_2)\Delta L,
\end{align*}
$$

where $\Delta L = L - L^*$ and<fnref n="1" />

$$
\frac{\partial \mathcal{H}_K}{\partial \Delta \Upsilon} = -3\eta \left\{ (p^2 C_2 + p(p+1) C_3) \Delta \Upsilon + p C_2 \Delta L^* \right\}. \quad (C.2)
$$

<footnote n="1">¹$\Delta L^*$ is defined by Eq. (4.26).</footnote>
## C.2 Expression of the matrix $Q_6$

We give here the matrix $Q_6$ that appears in Eq. (4.37). We denote

$$ 
\begin{aligned}
r_j &= \sqrt{2D_j}, \quad c_3 = \cos \sigma_3, \quad s_3 = \sin \sigma_3, \quad s = \sin \frac{p\pi}{3}, \\
c &= \cos \frac{p\pi}{3}, \quad I = 3\eta p (C_3 + p(C_2 + C_3)),
\end{aligned} 
\quad (C.3) $$

and obtain $Q_6 =$

$$ 
\begin{pmatrix}
I c s r_1^2 & I s r_1 r_2 & c_3 I s r_1 r_3 & I s^2 r_1^2 - \nu_3 & 0 & I s_3 s r_1 r_3 \\
0 & 0 & 0 & 0 & -\nu_3 & 0 \\
I s_3 c r_1 r_3 & I s_3 r_2 r_3 & c_3 I s_3 r_3^2 & I s_3 s r_1 r_3 & 0 & I s_3^2 r_3^2 - \nu_3 \\
-I c^2 r_1^2 + \nu_3 & -I c r_1 r_2 & -c_3 I c r_1 r_3 & -I c s r_1^2 & 0 & -I s_3 c r_1 r_3 \\
-I c r_1 r_2 & -I r_2^2 + \nu_3 & -c_3 I r_2 r_3 & -I s r_1 r_2 & 0 & -I s_3 r_2 r_3 \\
-c_3 I c r_1 r_3 & -c_3 I r_2 r_3 & -c_3^2 I r_3^2 + \nu_3 & -c_3 I s r_1 r_3 & 0 & -c_3 I s_3 r_3^2
\end{pmatrix} . 
\quad (C.4) $$

## C.3 Expression of the matrix $Q_0 + Q_1$

We give here the complete expression of the matrix $Q_0 + Q_1$ of the linearized system (6.24). We have

$$ Q_0 = \begin{pmatrix} \mathcal{Z}_0 & 0_{5 \times 2} \\ 0_{2 \times 5} & \mathcal{M}_0 \end{pmatrix} \quad \text{and} \quad Q_1 = \begin{pmatrix} \mathcal{Z}_1 & 0_{5 \times 2} \\ 0_{2 \times 5} & \mathcal{M}_1 \end{pmatrix}, \quad (C.5) $$

where

$$ \mathcal{M}_0 = \frac{27}{8} i \eta \begin{pmatrix} \frac{m_2}{m_0} & -\frac{m_2}{m_0} e^{i \pi / 3} \\ -\frac{m_1}{m_0} e^{-i \pi / 3} & \frac{m_1}{m_0} \end{pmatrix}, \quad (C.6) $$

$$ \mathcal{M}_1 = -\frac{21}{2} \eta \operatorname{diag} \left\{ \frac{q_1}{Q_1} \frac{m_0}{m_1} \left(1 - \frac{5}{7} i Q_1\right), \frac{q_2}{Q_2} \frac{m_0}{m_2} \left(1 - \frac{5}{7} i Q_2\right) \right\}, \quad (C.7) $$

$$ \mathcal{Z}_0 = \begin{pmatrix} 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -3\gamma\eta & 0 \\ 0 & 0 & 9i\eta\gamma^{-1/4} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{pmatrix}, \quad (C.8) $$

$$ \mathcal{Z}_1 = \begin{pmatrix} -\mathfrak{a}_1 & 0 & 0 & 3\gamma^{-1}\delta^{-1}\mathfrak{a}_1 & 3\gamma^{-1}\mathfrak{a}_1 \\ 0 & -\mathfrak{a}_2 & 0 & -3\gamma\delta\mathfrak{a}_2 & 3\gamma^{-1}\mathfrak{a}_2 \\ 0 & 0 & 0 & -\gamma[\delta\mathfrak{c}_2 + (1-\delta)\mathfrak{c}_1] & \gamma^{-1}(\mathfrak{c}_2 - \mathfrak{c}_1) \\ -(1-\delta)\mathfrak{b}_1 & \delta\mathfrak{b}_2 & 0 & 3\gamma[\delta^2\mathfrak{b}_2 + (1-\delta)^2\mathfrak{b}_1] & 3\gamma^{-1}[(1-\delta)\mathfrak{b}_1 - \delta\mathfrak{b}_2] \\ -\mathfrak{b}_1 & -\mathfrak{b}_2 & 0 & 3\gamma^{-1}\delta^{-1}\mathfrak{b}_1 - 3\gamma\delta\mathfrak{b}_2 & 3\gamma^{-1}(\mathfrak{b}_1 + \mathfrak{b}_2) \end{pmatrix}, \quad (C.9) $$

with

$$ 
\begin{aligned}
\delta &= \frac{m_1}{m_1 + m_2}, & \gamma &= \frac{m_1 + m_2}{m}, & \iota &= \frac{m_1 + m_2}{m_0} \\
\mathfrak{a}_j &= \frac{3\eta q_j}{\alpha_j Q_j} \gamma_j^{-2} \frac{m_0}{m_j}, & \mathfrak{b}_j &= 3\eta \frac{q_j}{Q_j} \frac{m_0}{m}, & \mathfrak{c}_j &= 78\eta q_j \frac{m_0}{m_j}.
\end{aligned} 
\quad (C.10) $$
Near $L_{4,5}$, the eigenvectors of $M_0 + M_1$, computed using results from Sect. 2.4.2, reveal that the Lagrange configuration corresponds to

$$ 
\begin{aligned}
\bar{\omega}_1 - \bar{\omega}_2 &= \frac{\pi}{3} + \frac{28 m_0^2 (m_1 q_2 / Q_2 + m_2 q_1 / Q_1)}{9 m_1 m_2 (m_1 + m_2)}, \\
\frac{e_1}{e_2} &= 1 + \frac{20 m_0^2 (q_2 m_1 - q_1 m_2)}{9 m_1 m_2 (m_1 + m_2)},
\end{aligned}
\qquad (C.11) $$

while the anti-Lagrange configuration complies with

$$ 
\begin{aligned}
\bar{\omega}_1 - \bar{\omega}_2 &= \frac{4\pi}{3} - \frac{28 m_0^2 (m_1 q_2 / Q_2 + m_2 q_1 / Q_1)}{9 m_1 m_2 (m_1 + m_2)}, \\
\frac{e_1}{e_2} &= \frac{m_2}{m_1} \left( 1 - \frac{20 m_0^2 (q_2 m_1 - q_1 m_2)}{9 m_1 m_2 (m_1 + m_2)} \right).
\end{aligned}
\qquad (C.12) $$

These last two equations are a perturbation due to tides of Eqs (3.44) & (3.45).

## C.4 Complete model for tidal interactions

The complete equations of motion governing the tidal evolution of a planar ($N+1$)-body system in a heliocentric reference frame, using a linear constant time-lag tidal model, are given, for $1 \le j \le N$, by (Mignard, 1979)

$$ 
\begin{aligned}
\ddot{\mathbf{r}}_j &= -\frac{\mu_j}{r_j^3} \mathbf{r}_j + \sum_{k \ne j} \mathcal{G} m_k \left( \frac{\mathbf{r}_k - \mathbf{r}_j}{|\mathbf{r}_k - \mathbf{r}_j|^3} - \frac{\mathbf{r}_k}{r_k^3} \right) + \frac{\mathbf{f}_j}{\beta_j} + \sum_{k \ne j} \frac{\mathbf{f}_k}{m_0}, \\
\ddot{\theta}_j &= -3 \Delta t_j \frac{\kappa_2^{(j)} \mathcal{G} m_0^2 R_j^3}{\alpha_j m_j r_j^8} \left[ \dot{\theta} r_j^2 - (\mathbf{r}_j \times \dot{\mathbf{r}}_j) \cdot \mathbf{k} \right],
\end{aligned}
\qquad (C.13) $$

where $\mathbf{r}_j$ is the heliocentric position vector, $\theta_j$ the rotation angle of the planet $j$, $\mathbf{k}$ is the unit vector normal to the orbital plane, and $\mathbf{f}_j$ is the force arising from the tidal potential energy created by the deformation of planet $j$ (Eq. (6.3))

$$ \mathbf{f}_j = -3 \frac{\kappa_2^{(j)} \mathcal{G} m_0^2 R_j^5}{r_j^8} \mathbf{r}_j - 3 \frac{\kappa_2^{(j)} \mathcal{G} m_0^2 R_j^5}{r_j^{10}} \Delta t_j \left[ 2 (\mathbf{r}_j \cdot \dot{\mathbf{r}}_j) \mathbf{r}_j + r_j^2 (\dot{\theta} \mathbf{r}_j \times \mathbf{k} + \dot{\mathbf{r}}_j) \right]. \qquad (C.14) $$
Bibliography

Adams, F. C. and Bloch, A. M. (2015). On the Stability of Extrasolar Planetary Systems and Other Closely Orbiting Pairs. *Monthly Notices of the Royal Astronomical Society*, 446, pp. 3676–3686.

Andoyer, H. (1923). *Cours de Mecanique Celeste*. Gauthier-Villars et Cie.

Barker, A. J. (2020). Tidal Dissipation in Evolving Low-Mass and Solar-Type Stars with Predictions for Planetary Orbital Decay. *Monthly Notices of the Royal Astronomical Society*, 498, pp. 2270–2294.

Beaugé, C., Michtchenko, T. A. and Ferraz-Mello, S. (2006). Planetary Migration and Extrasolar Planets in the 2/1 Mean-Motion Resonance. *Monthly Notices of the Royal Astronomical Society*, 365, pp. 1160–1170.

Beaugé, C. et al. (2007). Co-Orbital Terrestrial Planets in Exoplanetary Systems: A Formation Scenario. *Astronomy and Astrophysics*, Volume 463, Issue 1, February III 2007, pp.359-367, 463.1, p. 359.

Boué, G., Correia, A. C. M. and Laskar, J. (2019). On Tidal Theories and the Rotation of Viscous Bodies. 82, pp. 91–98.

Cardano, G. (1545). *Ars Magna*.

Christensen, R. M. (1982). *Theory of Viscoelasticity : An Introduction*. Second edition, United Kingdom edition. Academic press.

Claret, A. and Cunha, N. C. S. (1997). Circularization and Synchronization Times in Main-Sequence of Detached Eclipsing Binaries II. Using the Formalisms by Zahn. *Astronomy and Astrophysics*, v.318, p.187-197, 318, p. 187.

Correia, A. C. M., Bourrier, V. and Delisle, J.-B. (2020). Why Do Warm Neptunes Present Nonzero Eccentricity? *Astronomy & Astrophysics*, Volume 635, id.A37, 635, A37.

Correia, A. C. M. and Laskar, J. (2010). *Tidal Evolution of Exoplanets*. University of Arizona Press, Tucson, pp. 239–266.

Correia, A. C. M. (2009). Secular Evolution of a Satellite by Tidal Effect: Application to Triton. *The Astrophysical Journal*, 704, pp. L1–L4.

Correia, A. C. M. and Laskar, J. (2003). Different Tidal Torques on a Planet with a Dense Atmosphere and Consequences to the Spin Dynamics. *Journal of Geophysical Research (Planets)*, 108, p. 5123.

Correia, A. C. M. and Laskar, J. (2004). Mercury's Capture into the 3/2 Spin-Orbit Resonance as a Result of Its Chaotic Dynamics. *Nature*, 429, pp. 848–850.

Correia, A. C. M. et al. (2014). Deformation and Tidal Evolution of Close-in Planets and Satellites Using a Maxwell Viscoelastic Rheology. *Astronomy & Astrophysics*, Volume 571, id.A50, 571, A50.

Couetdic, J. et al. (2010). Dynamical Stability Analysis of the HD 202206 System and Constraints to the Planetary Orbits. *Astronomy and Astrophysics*, Volume 519, id.A10, 519, A10.
Couturier, J., Robutel, P. and Correia, A. C. M. (2022). Dynamics of Co-Orbital Exoplanets in a First-Order Resonance Chain with Tidal Dissipation. *Astronomy & Astrophysics*, 664, A1.

Couturier, J., Robutel, P. and Correia, A. C. M. (2021). An Analytical Model for Tidal Evolution in Co-Orbital Systems. I. Application to Exoplanets. *Celestial Mechanics and Dynamical Astronomy*, 133.8, p. 37.

Cresswell, P. and Nelson, R. P. (2008). Three-Dimensional Simulations of Multiple Protoplanets Embedded in a Protostellar Disc. *Astronomy and Astrophysics*, Volume 482, Issue 2, 2008, pp.677-690, 482.2, p. 677.

Darwin, G. H. (1880). On the Secular Changes in the Elements of the Orbit of a Satellite Revolving about a Tidally Distorted Planet. *Philosophical Transactions of the Royal Society of London Series I*, 171, pp. 713-891.

Delaunay, M. (1861). Théorie Du Mouvement de La Lune. *Monthly Notices of the Royal Astronomical Society*, 21, p. 80.

Delisle, J.-B. (2017). Analytical Model of Multi-Planetary Resonant Chains and Constraints on Migration Scenarios. *Astronomy & Astrophysics*, Volume 605, id.A96, 605, A96.

Delisle, J.-B. and Laskar, J. (2014). Tidal Dissipation and the Formation of Kepler Near-Resonant Planets. *Astronomy & Astrophysics*, Volume 570, id.L7, 570, p. L7.

Delisle, J.-B., Laskar, J. and Correia, A. C. M. (2014). Resonance Breaking Due to Dissipation in Planar Planetary Systems. *Astronomy & Astrophysics*, Volume 566, id.A137, 566, A137.

Dobrovolskis, A. R. and Lissauer, J. J. (2022). Do Tides Destabilize Trojan Exoplanets? *Icarus*, 385, p. 115087.

Efroimsky, M. (2012). Bodily Tides near Spin-Orbit Resonances. *Celestial Mechanics and Dynamical Astronomy*, 112, pp. 283-330.

Euler (1764). Considerationes de Motu Corporum Coelestium. Novi Commentarii Academiae Scientiarum Petropolitanae. *Berlin Acad*, 10, pp. 544-558.

Ford, E. B. and Gaudí, B. S. (2006). Observational Constraints on Trojans of Transiting Extrasolar Planets. *The Astrophysical Journal*, 652, pp. L137-L140.

Fortney, J. J. et al. (2018). *The Interior of Saturn*. Cambridge University Press, pp. 44–68.

Gascheau (1843). Examen d'une Classe d'équations Différentielles et Application à Un Cas Particulier Du Problème Des Trois Corps. *C. R. Acad. Sci. Paris*, 16.

Gastineau, M. and Laskar, J. (2011). TRIP: A Computer Algebra System Dedicated to Celestial Mechanics and Perturbation Series. *ACM Communications in Computer Algebra*, 44, pp. 194–197.

Giuppone, C. A., Benítez-Llambay, P. and Beaugé, C. (2012). Origin and Detectability of Co-Orbital Planets from Radial Velocity Data. *Monthly Notices of the Royal Astronomical Society*, 421, pp. 356–368.

Giuppone, C. A. et al. (2010). Dynamics of Two Planets in Co-Orbital Motion. *Monthly Notices of the Royal Astronomical Society*, 407, pp. 390–398.

Hara, N. C. et al. (2020). The SOPHIE Search for Northern Extrasolar Planets. XVI. HD 158259: A Compact Planetary System in a near-3:2 Mean Motion Resonance Chain. *Astronomy and Astrophysics*, 636, p. L6.

Henrard, J. and Lemaitre, A. (1983). A Second Fundamental Model for Resonance. *Celestial Mechanics*, 30, pp. 197–218.

Hill, G. W. (1913). Motion of a System of Material Points under the Action of Gravitation. *The Astronomical Journal*, 27, pp. 171–182.

Hut, P. (1980). Stability of Tidal Equilibrium. *Astronomy and Astrophysics*, vol. 92, no. 1-2, Dec. 1980, p. 167-170.

Lagrange (1772). Œuvres Complètes. Gauthier-Villars,
Lainey, V. (2016). Quantification of Tidal Parameters from Solar System Data. *Celestial Mechanics and Dynamical Astronomy*, 126, pp. 145–156.

Lainey, V. et al. (2012). Strong Tidal Dissipation in Saturn and Constraints on Enceladus’ Thermal State from Astrometry. *The Astrophysical Journal*, 752, p. 14.

Landau and Lifshitz (1980). *The Classical Theory of Fields*. Fourth. Vol. 2. Butterworth-Heinemann.

Laskar, J. (1990). The Chaotic Motion of the Solar System: A Numerical Estimate of the Size of the Chaotic Zones. *Icarus*, 88, pp. 266–291.

Laskar, J. (2000). On the Spacing of Planetary Systems. *Physical Review Letters*, 84, pp. 3240–3243.

Laskar, J., Boué, G. and Correia, A. C. M. (2012). Tidal Dissipation in Multi-Planet Systems and Constraints on Orbit Fitting. *Astronomy and Astrophysics*, 538, A105.

Laskar, J. (1993). Frequency Analysis for Multi-Dimensional Systems. Global Dynamics and Diffusion. *Physica D Nonlinear Phenomena*, 67, pp. 257–281.

Laskar, J. (2005). Note on the Generalized Hansen and Laplace Coefficients. *Celestial Mechanics and Dynamical Astronomy*, 91, pp. 351–356.

Laskar, J. (2017). Andoyer Construction for Hill and Delaunay Variables. *Celestial Mechanics and Dynamical Astronomy*, 128, pp. 475–482.

Laskar, J. and Robutel, P. (1995). Stability of the Planetary Three-Body Problem. I. Expansion of the Planetary Hamiltonian. *Celestial Mechanics and Dynamical Astronomy*, 62, pp. 193–217.

Laughlin, G. and Chambers, J. E. (2002). Extrasolar Trojans: The Viability and Detectability of Planets in the 1:1 Resonance. *The Astronomical Journal*, 124, pp. 592–600.

Leleu, A., Robutel, P. and Correia, A. C. M. (2018). On the Coplanar Eccentric Non-Restricted Co-Orbital Dynamics. *Celestial Mechanics and Dynamical Astronomy*, 130, p. 24.

Leleu, A. et al. (2017). Detection of Co-Orbital Planets by Combining Transit and Radial-Velocity Measurements. *Astronomy & Astrophysics*, Volume 599, id.L7, 599, p. L7.

Leleu, A. et al. (2021). Alleviating the Transit Timing Variation Bias in Transit Surveys. I. RIVERS: Method and Detection of a Pair of Resonant Super-Earths around Kepler-1705. *Astronomy & Astrophysics*, Volume 655, id.A66, 655, A66.

Leleu, A. et al. (2022). Alleviating the Transit Timing Variation Bias in Transit Surveys. II. RIVERS: Twin Resonant Earth-sized Planets around Kepler-1972 Recovered from a Kepler False Positive. *Astronomy & Astrophysics*, Volume 661, id.A141, 661, A141.

Leleu, A., Coleman, G. A. L. and Ataiee, S. (2019). Stability of the Co-Orbital Resonance under Dissipation. Application to Its Evolution in Protoplanetary Discs. *Astronomy and Astrophysics*, 631, A6.

Leleu, A., Robutel, P. and Correia, A. C. M. (2015). Detectability of Quasi-Circular Co-Orbital Planets. Application to the Radial Velocity Technique. *Astronomy and Astrophysics*, 581, A128.

Lillo-Box, J. et al. (2018). The TROY Project: Searching for Co-Orbital Bodies to Known Planets. I. Project Goals and First Results from Archival Radial Velocity. *Astronomy & Astrophysics*, Volume 609, id.A96, 609, A96.

Luhn, J. K. et al. (2019). Retired A Stars and Their Companions. VIII. 15 New Planetary Signals around Subgiants and Transit Parameters for California Planet Search Planets with Subgiant Hosts. *The Astronomical Journal*, 157, p. 149.

Lyra, W. et al. (2009). Standing on the Shoulders of Giants. Trojan Earths and Vortex Trapping in Low Mass Self-Gravitating Protoplanetary Disks of Gas and Solids. *Astronomy and Astrophysics*, Volume 493, Issue 3, 2009, pp.1125-1139, 493.3, p. 1125.

MacDonald, G. J. F. (1964). Tidal Friction. *Reviews of Geophysics and Space Physics*, 2, pp. 467–541.
Malhotra, R. (2022). New Results on Orbital Resonances. *IAU Symposium*, 364, pp. 85–101.

Meyer, K. R. and Hall, G. R. (2009). *Introduction to Hamiltonian Dynamical Systems and the N-Body Problem*. Second edition. Vol. 90. Springer.

Michtchenko, T. A., Beaugé, C. and Ferraz-Mello, S. (2008). Dynamic Portrait of the Planetary 2/1 Mean-Motion Resonance - I. Systems with a More Massive Outer Planet. *Monthly Notices of the Royal Astronomical Society*, 387, pp. 747–758.

Mignard, F. (1979). The Evolution of the Lunar Orbit Revisited. I. *Moon and Planets*, 20, pp. 301–315.

Moeckel, R. (2017). Minimal Energy Configurations of Gravitationally Interacting Rigid Bodies. *Celestial Mechanics and Dynamical Astronomy*, 128, pp. 3–18.

Morbidelli, A. (2002). *Modern Celestial Mechanics: Aspects of Solar System Dynamics*.

Namouni, F. (1999). Secular Interactions of Coorbiting Objects. *Icarus*, 137, pp. 293–314.

Niederman, L., Pousse, A. and Robutel, P. (2020). On the Co-orbital Motion in the Three-Body Problem: Existence of Quasi-periodic Horseshoe-Shaped Orbits. *Communications in Mathematical Physics*, 377, pp. 551–612.

Nimmo, F. and Schenk, P. (2006). Normal Faulting on Europa: Implications for Ice Shell Properties. *Journal of Structural Geology*, 28, pp. 2194–2203.

Petit, A. C. et al. (2020). The Path to Instability in Compact Multi-Planetary Systems. *Astronomy & Astrophysics*, Volume 641, id.A176, 641, A176.

Pierens, A. and Raymond, S. N. (2014). Disruption of Co-Orbital (1:1) Planetary Resonances during Gas-Driven Orbital Migration. *Monthly Notices of the Royal Astronomical Society*, 442, pp. 2296–2303.

Pousse, A., Robutel, P. and Vienne, A. (2017). On the Co-Orbital Motion in the Planar Restricted Three-Body Problem: The Quasi-Satellite Motion Revisited. *Celestial Mechanics and Dynamical Astronomy*, 128, pp. 383–407.

Pucacco, G. (2021). Normal Forms for the Laplace Resonance. *Celestial Mechanics and Dynamical Astronomy*, 133, p. 11.

Quillen, A. C. et al. (2017). Obliquity Evolution of the Minor Satellites of Pluto and Charon. *Icarus*, 293, pp. 94–113.

Ragazzo, C. and Ruiz, L. S. (2015). Dynamics of an Isolated, Viscoelastic, Self-Gravitating Body. *Celestial Mechanics and Dynamical Astronomy*, 122, pp. 303–332.

Roberts, J. H. and Nimmo, F. (2008). Tidal Heating and the Long-Term Stability of a Subsurface Ocean on Enceladus. *Icarus*, 194, pp. 675–689.

Robutel, P. and Gabern, F. (2006). The Resonant Structure of Jupiter's Trojan Asteroids - I. Long-term Stability and Diffusion. *Monthly Notices of the Royal Astronomical Society*, 372, pp. 1463–1482.

Robutel, P. and Pousse, A. (2013). On the Co-Orbital Motion of Two Planets in Quasi-Circular Orbits. *Celestial Mechanics and Dynamical Astronomy*, 117, pp. 17–40.

Rodríguez, A., Giuppone, C. A. and Michtchenko, T. A. (2013). Tidal Evolution of Close-in Exoplanets in Co-Orbital Configurations. *Celestial Mechanics and Dynamical Astronomy*, 117, pp. 59–74.

Segatz, M. et al. (1988). Tidal Dissipation, Surface Heat Flow, and Figure of Viscoelastic Models of Io. *Icarus*, 75, pp. 187–206.

Vokrouhlický, D. and Nesvorný, D. (2014). Transit Timing Variations for Planets Co-orbiting in the Horseshoe Regime. *The Astrophysical Journal*, 791, p. 6.

Wahr, J. (1996). *Geodesy and Gravity*. Samizdat Press.

Yoder, C. F. et al. (1983). Theory of Motion of Saturn's Coorbiting Satellites. *Icarus*, 53, pp. 431–443.
<figure><img src="image_31.png" /></figure>
# RÉSUMÉ

Les systèmes planétaires peuvent adopter des configurations remarquables. L'une d'elles, dite co-orbitale, se produit lorsque deux planètes ont la même période orbitale autour de leur étoile, c'est à dire, quand elles sont en résonance de moyen mouvement 1 : 1. Même au sein de la résonance co-orbitale, de nombreuses trajectoires sont possibles. Les plus simples sont connues depuis le 18ème siècle, comme deux corps co-orbitaux sur des orbites planes et circulaires formant avec leur étoile un triangle équilatéral tournant. Cependant, certaines configurations plus subtiles n'ont été découvertes que récemment. Dans le cas de deux points matériels à faibles excentricités et inclinaisons, les aspects analytiques de la dynamique sont bien compris. Cependant, de grandes excentricités ou inclinaisons sont responsables de changements topologiques dans l'espace des phases, tandis que des corps étendus peuvent engendrer la dissipation de l'énergie mécanique, et la dynamique dans ces cas a encore des zones d'ombre.

Aucune des huit planètes du système Solaire ne co-orbitent ensemble le Soleil, bien que des corps co-orbitaux existent dans le système Solaire, soit entre deux objets mineurs (orbitant une planète), soit entre une planète et un objet mineur (orbitant le Soleil). Cette absence de planètes co-orbitales n'est a priori pas la norme dans les systèmes exoplanétaires, puisque les modèles de formation prédisent leur existence. Pourtant, en dépit de milliers de détections d'exoplanètes, aucune paire de planètes co-orbitales n'a été détectée à ce jour. Bien que cela puisse s'expliquer par des biais observationnels, nous montrons dans ce manuscrit que les forces de marées sont responsables de la destruction des paires de planètes co-orbitales. Nous construisons un modèle analytique de marées du système plan étoile-planète-planète, basé sur une extension, avec dissipation de marées, du formalisme Hamiltonien. Le modèle fournit une expression analytique précise de la durée de vie de la paire, dépendant des paramètres, et qui permet de prédire quelles exoplanètes déjà découvertes pourraient avoir un compagnon co-orbital non détecté.

Les modèles de formation prédisent aussi qu'un nombre important de planètes co-orbitales sont formées au sein d'une chaîne de résonance. Ainsi, nous étendons l'étude précédente au cas où la paire est au sein d'une chaîne de résonance. Plus précisement, nous construisons un modèle Hamiltonien de la chaîne de résonance *p* : *p* : *p* + 1 où la paire de co-orbitaux est en résonance de moyen mouvement du premier ordre *p* : *p* + 1 avec une troisième planète externe, *p* étant un petit entier. Après comparaison des familles d'équilibres du modèle avec les familles d'orbites quasipériodiques correspondantes dans le système complet, nous ajoutons la dissipation de marées au modèle à l'aide d'un formalisme pseudo-Hamiltonien. Nous montrons que cette chaîne de résonance met en scène une résonance séculaire 1 : 1 entre la libration de l'angle co-orbital et la précession des péricentres, et en analysant les valeurs propres du système différentiel linéarisé au voisinage des familles d'équilibres, nous montrons comment les marées stabilisent le système aux alentours de cette résonance séculaire, rendant la paire de planètes co-orbitales bien plus stable quand elle est dans la chaîne de résonance *p* : *p* : *p* + 1.

# MOTS CLÉS

Mécanique Céleste ○ Marées ○ Problème des Trois & Quatre Corps ○ Résonance ○ Co-orbital

# ABSTRACT

Planetary systems can adopt remarkable configurations. One of them, said co-orbital, occurs when two planets share the same orbital period around their star, that is, when they are in a 1 : 1 mean motion resonance. Even within the co-orbital motion, many trajectories are possible. The simplest ones are known since the 18th century, like two co-orbital bodies on planar and circular orbits forming with the star a rotating equilateral triangle. However, some more subtle configurations were not discovered until recently. In the case of point mass bodies with small eccentricities and inclinations, the analytical features of the dynamics are well understood. However, large eccentricities and inclinations are responsible for topological changes in the phase space, while extended bodies can lead to dissipation of the orbital energy, and the dynamics in these cases still has grey areas.

None of the eight planets of the Solar system are in co-orbital motion together, although co-orbital bodies exist in the Solar system, either between two minor objects (orbiting a planet) or between a planet and a minor object (orbiting the Sun). This absence of planet-planet co-orbital motion should not be the norm in exoplanetary systems, as the formation models predict their existence. Nevertheless, despite thousands of exoplanets detection, no pair of co-orbital planets has been detected so far. While this may be in part explained by detection bias, we show in this manuscript that tidal effects are responsible for the disruption of co-orbital pairs of planets. We build an analytical tidal model of the planar system star-planet-planet, based on an extension of the point-mass Hamiltonian formalism with tidal dissipation. The model provides an accurate analytical expression for the lifetime of the pair depending on the parameters, allowing to predict which already discovered exoplanets might have an undetected co-orbital companion.

Formation models also predict that a significant number of co-orbital planets are formed within a resonance chain. Therefore, we extend the previous study to the case where the pair is within a resonance chain. More precisely, we build a Hamiltonian model of the resonance chain *p* : *p* : *p* + 1 where the co-orbital pair is in a first-order mean motion resonance *p* : *p* + 1 with an outermost third planet, *p* being a small integer. After comparing the families of equilibria of the model with the associated families of periodic orbits of the complete system, we add tidal dissipation to the model using a pseudo-Hamiltonian formalism. We show that this resonance chain features a 1 : 1 secular resonance between the libration of the co-orbital angle and the precession of the pericentres, and by analizing the eigenvalues of the linearized system in the vicinity of the families of equilibria, we show how tides stabilize the system around this secular resonance, making co-orbital pairs of exoplanets much more stable while inside the *p* : *p* : *p* + 1 resonance chain.

# KEYWORDS