# Dynamics of co-orbital exoplanets in a first-order resonance chain with tidal dissipation

Jérémy Couturier, Philippe Robutel, Alexandre C.M. Correia

► To cite this version:

Jérémy Couturier, Philippe Robutel, Alexandre C.M. Correia. Dynamics of co-orbital exoplanets in a first-order resonance chain with tidal dissipation. Astronomy and Astrophysics - A&A, 2022, 664, pp.A1. 10.1051/0004-6361/202243261. [hal-03835319](https://HALIdHAL-03835319)

HAL Id: hal-03835319

https://cnrs.hal.science/hal-03835319

Submitted on 31 Oct 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Dynamics of co-orbital exoplanets in a first-order resonance chain with tidal dissipation

Jérémy Couturier<fnref n="1" />,*, Philippe Robutel<fnref n="1" />,*, and Alexandre C.M. Correia<fnref n="2" />,<fnref n="1" />,†

Received 3 February 2022 / Accepted 8 April 2022

## ABSTRACT

Co-orbital planets (in a 1 : 1 mean motion resonance) can be formed within a Laplace resonance chain. We develop a secular model to study the dynamics of the resonance chain *p*: *p*: *p* + 1, where the co-orbital pair is in a first-order mean motion resonance with the outermost third planet. Our model takes into account tidal dissipation through the use of a Hamiltonian version of the constant time-lag model, which extends the Hamiltonian formalism of the point-mass case. We show the existence of several families of equilibria, and how these equilibria extend to the complete system. In one family, which we call the main branch, a secular resonance between the libration frequency of the co-orbitals and the precession frequency of the pericentres has unexpected dynamical consequences when tidal dissipation is added. We report the existence of two distinct mechanisms that make co-orbital planets much more stable within the *p*: *p*: *p* + 1 resonance chain rather than outside it. The first is due to negative real parts of the eigenvalues of the linearised system with tides, in the region of the secular resonance mentioned above. The second comes from non-linear contributions of the vector field and is due to eccentricity damping. These two stabilising mechanisms increase the chances of a future detection of exoplanets in the co-orbital configuration.

**Key words.** celestial mechanics

## 1. Introduction

Co-orbital systems have been extensively studied since the discovery of equilibria in the three-body problem by Euler (1764) and Lagrange (1772). In the hierarchical case, that is when $(m_1 + m_2)/m_0 < 1/27$, where $m_0$ is the mass of the central body and $m_1$ and $m_2$ are the co-orbital masses, Gascheau (1843) showed that the equilateral equilibria, where the bodies are at the vertices of an equilateral triangle, are linearly stable. This result, combined with the discovery of several examples of co-orbital bodies, such as the Jovian trojans or the horseshoe-shaped orbits of Janus and Epimetheus around Saturn, contributed to increasing the interest of scientists for systems of this kind.

In the planar and circular case, for masses in the range $3 \times 10^{-4} < (m_1 + m_2)/m_0 < 1/27$, the angle $\lambda_1 - \lambda_2$ between the co-orbitals librates around its equilibrium of $\pm 60^\circ$ in trajectories commonly called tadpole orbits, where the co-orbital angle is bounded by<fnref n="1" /> $23.9^\circ < \lambda_1 - \lambda_2 < 180^\circ$. However, for masses $(m_1 + m_2)/m_0 < 3 \times 10^{-4}$ (e.g. Laughlin & Chambers 2002), a separatrix in the phase space delimits a region of stable trajectories of another kind, generally called horseshoe orbits, where the critical angle $\lambda_1 - \lambda_2$ librates around $180^\circ$ with at least $312.2^\circ$ of amplitude. For low eccentricities and small libration amplitudes, still in the planar case, the existence of two proper modes called Lagrange and anti-Lagrange have been shown, numerically by Giuppone et al. (2010) and analytically by Robutel & Pousse (2013). In the Lagrange configuration, the pericentres of the orbits do not precess and verify the relation $\bar{\omega}_1 - \bar{\omega}_2 = 60^\circ$,

whereas in the anti-Lagrange configuration, both orbits precess at the same frequency while maintaining the relation $\bar{\omega}_1 - \bar{\omega}_2 = 240^\circ$. For low eccentricities but very large libration amplitudes, in the region of horseshoe-shaped orbits, Couturier et al. (2021) showed that the Lagrange and anti-Lagrange configurations correspond to $\bar{\omega}_1 - \bar{\omega}_2 = 0^\circ$ and $\bar{\omega}_1 - \bar{\omega}_2 = 180^\circ$, respectively. High eccentricities give rise to topological changes in the phase space (see Leleu et al. 2018), and thus to many more exotic trajectories, while the dynamics of the inclined problem is even more complex by allowing, among other things, transitions between these orbits and the retrograde co-orbitals (Namouni 1999).

The discovery of exoplanets raised the question of the existence of co-orbital planets, which are absent from the Solar System. Accretion in situ at the equilateral equilibria of a primary or capture in the 1:1 resonance of planets formed in other parts of the system are two possible scenarios of formation of such systems (e.g. Laughlin & Chambers 2002; Cresswell & Nelson 2009). The stability of co-orbital planets formed in a disc has been studied by Leleu et al. (2019), who show that under dissipative interactions with the gas disc, the equilateral equilibria can be either attractive or repulsive, depending on the co-orbital mass ratio and the parameters of the disc. Moreover, Leleu et al. (2019) show that, at least around low-mass stars, co-orbital exoplanets generally end up in a tadpole configuration and often within a Laplace resonance chain.

For co-orbital exoplanets orbiting close to the host star, tidal dissipation induced by the differential gravitational interaction leads to a long-term evolution of the orbits. Couturier et al. (2021) shows that, for a pair of co-orbital exoplanets orbiting a star, the equilateral Lagrangian equilibria are always repulsive under tidal interactions, and that regardless of the parameters of

<footnote n="1">¹ IMCCE, UMR8028 CNRS, Observatoire de Paris, PSL Univ., Sorbonne Univ., 77 av. Denfert-Rochereau, 75014 Paris, France
e-mail: jeremy.couturier@obspm.fr</footnote>

<footnote n="2">² CFisUC, Departamento de Física, Universidade de Coimbra, 3004-516 Coimbra, Portugal</footnote>

<footnote n="1">¹ The exact value of the lower bound is $2 \arcsin[(\sqrt{2} - 1)/2] \approx 23.9^\circ$ (e.g. Robutel & Pousse 2013).</footnote>
the system, the destruction of the co-orbital motion is unavoidable. However, the discovery of co-orbital exoplanets is still possible because the destruction time is strongly dependent on the parameters of the system and can easily be greater than the lifetime of the host star. Couturier et al. (2021) neglected any interaction with possible other planets in the system. In this paper we extend this work to the case where the pair of co-orbital exoplanets interacts with an outermost third planet, in a first-order mean motion resonance with the co-orbitals. More precisely, we study the Laplace resonance chain $p : p : p + 1$, where $p$ is a small integer.

In Sect. 2, we study the point-mass $p : p : p + 1$ resonance chain in the absence of tides. We show how rich and complex the dynamics of this chain is, and we conclude the section with the presentation of the stability map of the chain. In Sect. 3, we include tidal dissipation in the model by an extension of the Hamiltonian formalism. We study the linearised system in the vicinity of the equilibria, and by computing the real parts of the eigenvalues, we show the existence of a linearly stable zone around the 1:1 secular resonance between the libration frequency of the co-orbital and the precession frequency of the pericentres. In Sect. 4, we compare the analytical results with numerical simulations. They confirm the results of Sects. 2 and 3 and highlight the existence of a stabilisation mechanism of the co-orbitals due to eccentricity damping. We discuss our results in Sect. 5. In Table A.1, we list the notations used throughout this paper. Appendix F completes Sect. 4 with more numerical simulations and a complete discussion on the influence of the mass of the third planet on the co-orbital dynamics.

## 2. The $p : p : p + 1$ resonance chain

### 2.1. The Hamiltonian of the problem

In this section, we study an occurrence of the point-mass planar four-body problem. We construct the Hamiltonian associated with the resonance chain $p : p : p + 1$, where a central body, a star of mass $m_0$, is orbited by two co-orbital planets of mass $m_1$ and $m_2$ and a third planet of mass $m_3$ further away from the star and in a first-order mean motion resonance with the pair of co-orbital planets. Although we write all equations for a general value of the integer $p$, the figures are restricted to the case $p = 1$ where the nominal period of the third planet is twice that of the co-orbitals. For all planets we define the quantities $\beta_j = m_0 m_j / (m_0 + m_j)$ and $\mu_j = G(m_0 + m_j)$, where $G$ is the gravitational constant.

#### 2.1.1. The averaged Hamiltonian

In order to define a canonical coordinate system related to the semi-major axis $a_j$, the eccentricity $e_j$, the mean longitude $\lambda_j$, and the longitude of the pericentre $\varpi_j$ of planet $j$, we first consider Poincaré heliocentric coordinates $(\tilde{\Lambda}_j, \lambda_j, \tilde{D}_j, -\varpi_j)$, where

$$ \tilde{\Lambda}_j = \beta_j \sqrt{\mu_j a_j} \quad \text{and} \quad \tilde{D}_j = \tilde{\Lambda}_j \left(1 - \sqrt{1 - e_j^2}\right). \qquad (1) $$

In these coordinates, the Hamiltonian derives from the symplectic form

$$ \Omega = \sum_{j \in \{1, 2, 3\}} (d\lambda_j \wedge d\tilde{\Lambda}_j - d\varpi_j \wedge d\tilde{D}_j). \qquad (2) $$

Following Laskar & Robutel (1995), the planetary Hamiltonian is written

$$ H = H_K(\tilde{\Lambda}_j) + H_P(\tilde{\Lambda}_j, \lambda_j, \tilde{D}_j, \varpi_j), \qquad (3) $$

where the Keplerian part, due to star-planet interactions, reads

$$ H_K = - \sum_{j \in \{1, 2, 3\}} \frac{\beta_j^3 \mu_j^2}{2 \tilde{\Lambda}_j^2}, \qquad (4) $$

and the perturbation $H_P$, whose size relative to $H_K$ is of order $(m_1 + m_2 + m_3)/m_0$, is expanded in power series of the eccentricities. We assume that the system is close to the resonance $p : p : p + 1$. This means that the nominal mean motions verify

$$ n_{1,0} = n_{2,0} = \eta = \frac{p+1}{p} n_{3,0}, \qquad (5) $$

while the nominal semi-major axes $a_{j,0}$ are related to $n_{j,0}$ by the Kepler law $n_{j,0}^2 a_{j,0}^3 = \mu_0 = Gm_0$. The $a_j$ are always close to their nominal value$^2$ $a_{j,0}$, and the $\tilde{\Lambda}_j$ stay close to the quantities $\Lambda_j^*$ defined as

$$ \Lambda_j^* = m_j \sqrt{\mu_0 a_{j,0}}. \qquad (6) $$

To study the dynamics in the vicinity of the resonance, we expand the Hamiltonian in the neighbourhood of $(\Lambda_1^*, \Lambda_2^*, \Lambda_3^*)$. An expansion at order 2 in the Keplerian part and at order 0 in the perturbative part generates remainders of the same size, and we limit ourselves to

$$ H_P(\tilde{\Lambda}_j, \lambda_j, \tilde{D}_j, \varpi_j) = H_P(\Lambda_j^*, \lambda_j, \tilde{D}_j, \varpi_j). \qquad (7) $$

A suitable linear change of variables to deal with the $p : p : p + 1$ resonance chain is (e.g. Delisle 2017)

$$ \begin{pmatrix} \xi \\ \xi_2 \\ \xi_3 \\ \sigma_1 \\ \sigma_2 \\ \sigma_3 \end{pmatrix} = \begin{pmatrix} 1 & -1 & 0 & 0 & 0 & 0 \\ 0 & p & -p & 0 & 0 & 0 \\ 0 & -p & p+1 & 0 & 0 & 0 \\ 0 & -p & p+1 & 1 & 0 & 0 \\ 0 & -p & p+1 & 0 & 1 & 0 \\ 0 & -p & p+1 & 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} \lambda_1 \\ \lambda_2 \\ \lambda_3 \\ -\varpi_1 \\ -\varpi_2 \\ -\varpi_3 \end{pmatrix}, \qquad (8) $$

which is canonical if we transform the actions according to

$$ \begin{align} (\tilde{\Lambda}_1, \tilde{\Lambda}_2, \tilde{\Lambda}_3, \tilde{D}_j) &\mapsto (L', \Gamma', G', D'_j) \nonumber \\ &= \left( \tilde{\Lambda}_1, \frac{p+1}{p} (\tilde{\Lambda}_1 + \tilde{\Lambda}_2) + \tilde{\Lambda}_3, \sum_{j \le 3} (\tilde{\Lambda}_j - \tilde{D}_j), \tilde{D}_j \right). \tag{9} \end{align} $$

Since the total angular momentum $G'$ is a first integral, the Hamiltonian does not depend on the angle $\xi_3$. Moreover, in the $p : p : p + 1$ resonance, the angle $\xi_2$ is fast circulating and we average over it. The averaged Hamiltonian reads

$$ \begin{align} H'(L', \Gamma', G', \xi, D'_j, \sigma_j) \nonumber \\ &= H'_K(L', \Gamma', G', D'_j) + \frac{1}{2\pi} \int_0^{2\pi} H'_P(\xi, \xi_2, D'_j, \sigma_j) d\xi_2. \tag{10} \end{align} $$

<footnote>$^2$ The nominal semi-major axes are defined with $\mu_0$ instead of $\mu_j$, which conveniently yields $a_{1,0} = a_{2,0} := \bar{a}$. This approximation is valid since the subsequent error is of order $O(m_j/m_0)$, while the width of the resonance is of order $O(\sqrt{m_j/m_0})$.</footnote>
This change of variable, along with the averaging process, allows the two degrees of freedom associated with ($\xi_2, \xi_3, \Gamma', G'$) to be lost, and we are left with four degrees of freedom. After the averaging process, the scaling factor $\Gamma'$ and the angular momentum $G'$ are both parameters and a rescaling by $\Gamma'$ reduces the dependency to only one parameter. As we study the effect of tidal dissipation on the dynamics in Sect. 3, it is actually more convenient to normalise by the constant $\Gamma^* = (p+1)(\Lambda_1^* + \Lambda_2^*)/p + \Lambda_3^*$ rather than by $\Gamma'$, which is not constant when dissipation is present. In other words, we perform the canonical transformation

$$ \mathcal{H} = \frac{H'}{\Gamma^*}, \quad L = \frac{L'}{\Gamma^*}, \quad G = \frac{G'}{\Gamma^*}, \quad \Gamma = \frac{\Gamma'}{\Gamma^*}, \quad D_j = \frac{D'_j}{\Gamma^*}, \quad \Lambda_j = \frac{\tilde{\Lambda}_j}{\Gamma^*}, \quad (11) $$

while the angles are unchanged.

### 2.1.2. Expansion of the Keplerian part

As mentioned above, the Keplerian part of the Hamiltonian is expanded at second order in the vicinity of the $\Lambda_j^*$. If we denote $\Delta\tilde{\Lambda}_j = \tilde{\Lambda}_j - \Lambda_j^*$, the expansion reads

$$ H_K = \sum_{j=1}^{3} n_{j,0} \Delta \tilde{\Lambda}_j - \frac{3}{2} \sum_{j=1}^{3} \frac{n_{j,0}}{\Lambda_j^*} \Delta \tilde{\Lambda}_j^2. \quad (12) $$

Substituting $\tilde{\Lambda}_j - \Lambda_j^*$ for $\Delta\tilde{\Lambda}_j$, we obtain

$$ H_K = 4 \sum_{j=1}^{3} n_{j,0} \tilde{\Lambda}_j - \frac{3}{2} \sum_{j=1}^{3} \frac{n_{j,0}}{\Lambda_j^*} \tilde{\Lambda}_j^2 - \frac{5}{2} \sum_{j=1}^{3} n_{j,0} \Lambda_j^*. \quad (13) $$

The third term is constant and can be removed without changing the dynamics. Performing the normalisation (Eq. (11)) and the change of variable (Eq. (9)), we have

$$ \mathcal{H}_K = -\frac{3}{2}\eta \left\{ C_1 L^2 + C_2 (p(\Gamma - \Upsilon) - L)^2 + C_3 p(p+1) \left(\Upsilon - \frac{p\Gamma}{p+1}\right)^2 \right\} + \frac{4\eta p \Gamma}{p+1}, \quad (14) $$

where we denoted $\Upsilon = G + D_1 + D_2 + D_3 = \sum_j \Lambda_j$ and $C_j = \Gamma^*/\Lambda_j^*$, that is

$$ \begin{aligned} C_1 &= \left(\frac{p+1}{p}\right)^{1/3} \frac{m_3}{m_1} + \frac{p+1}{p} \left(1 + \frac{m_2}{m_1}\right), \\ C_2 &= \frac{m_1}{m_2} C_1, \quad C_3 = 1 + \left(\frac{p+1}{p}\right)^{2/3} \frac{m_1 + m_2}{m_3}. \end{aligned} \quad (15) $$

Without dissipation, $\Gamma$ is constant and almost equal to 1 and we simply evaluate $\mathcal{H}_K$ at $\Gamma = 1$, hence achieving the reduction to only one parameter<fnref n="3" />. In this case, the last term is constant and can also be removed.

Instead of the variables $L, G$, and $\Gamma$, we can use the variables $\Delta L, \Delta G$, and $\Delta \Gamma$ defined by their difference to the Keplerian

resonance (Eq. (6)). In this case, the approximation $\Gamma = 1$ becomes $\Delta\Gamma = 0$ and Eq. (12) yields, once normalised,

$$ \mathcal{H}_K = -\frac{3}{2}\eta \left( (p^2 C_2 + p(p+1)C_3) \Delta\Upsilon^2 + 2pC_2\Delta\Upsilon\Delta L + (C_1 + C_2)\Delta L^2 \right), \quad (16) $$

where $\Delta\Upsilon = \Delta G + D_1 + D_2 + D_3 = \sum_j \Delta\Lambda_j$. We find the Hamiltonian (16) to be well adapted to the analytical work derived in Sect. 2.2.2, while we use the Hamiltonian (14) in the remaining sections. Moreover, we do not perform the evaluation $\Gamma = 1$ in Sect. 3, where tidal dissipation is present and $\Gamma$ is a variable quantity. Both Hamiltonians yield the same dynamics and it is only a matter of preference.

### 2.1.3. Expansion of the perturbative part

The perturbative part $\mathcal{H}_P$ of the Hamiltonian is expanded in power series of the eccentricities. To this end, we separate the contributions due to the interactions between each pair of planets as

$$ \mathcal{H}_P = \mathcal{H}_{1,2} + \mathcal{H}_{1,3} + \mathcal{H}_{2,3}. \quad (17) $$

We denote $X_j = \sqrt{2\tilde{D}_j/\tilde{\Lambda}_j} e^{i\varpi_j} = e_j e^{i\varpi} + O(e_j^3)$ For a pair $(p_1, p_2) \in \{(1, 2), (1, 3), (2, 3)\}$ of planets, the perturbation to the Hamiltonian due to their mutual interaction reads (Laskar & Robutel 1995)

$$ \mathcal{H}_{p_1, p_2} = \sum_{\mathbf{k} \in \mathbb{Z}^2} \left( \sum_{\mathbf{q} \in \mathbb{N}^4} \Psi_k^q X_{p_1}^{q_1} X_{p_2}^{q_2} \bar{X}_{p_1}^{q_3} \bar{X}_{p_2}^{q_4} \right) e^{i(k_1 \lambda_{p_1} + k_2 \lambda_{p_2})}. \quad (18) $$

For a non-zero $\Psi_k^q$, the conservation of the angular momentum imposes on the tuples $\mathbf{q} = (q_1, q_2, q_3, q_4) \in \mathbb{N}^4$ and $\mathbf{k} = (k_1, k_2) \in \mathbb{Z}^2$ to verify the d'Alembert rule:

$$ k_1 + k_2 + q_1 + q_2 - q_3 - q_4 = 0. \quad (19) $$

This rule, combined with the averaging process, implies that $\mathcal{H}_{1,2}$ has no odd term in eccentricity, while $\mathcal{H}_{1,3}$ and $\mathcal{H}_{2,3}$ have no term of order 0. Since we limit ourselves to the second order in eccentricity, we write

$$ \mathcal{H}_{1,2} = \mathcal{H}^{(0)} + \mathcal{H}_{1,2}^{(2)}, \quad \mathcal{H}_{j,3} = \mathcal{H}_{j,3}^{(1)} + \mathcal{H}_{j,3}^{(2)}, \quad (20) $$

where the superscript refers to the order in eccentricity, while the subscript refers to the considered pair of planets. The Hamiltonian $\mathcal{H}^{(0)}$ has no subscript since only the pair of co-orbitals yields terms of order 0 and no confusion is possible.

Following Laskar & Robutel (1995), $\mathcal{H}_{1,3}$ and $\mathcal{H}_{2,3}$ can be written

$$ \mathcal{H}_{j,3}^{(1)} = \frac{m_j n_{3,0}}{m_0 C_3} \left\{ C_{p,1}^{(1)} \sqrt{2C_j D_j} \cos \left( p\delta_{j,1}\xi - \sigma_j \right) + C_{p,2}^{(1)} \sqrt{2C_j D_3} \cos \left( p\delta_{j,1}\xi - \sigma_3 \right) \right\} \quad (21) $$

and

$$ \begin{aligned} \mathcal{H}_{j,3}^{(2)} = \frac{2m_j n_{3,0}}{m_0 C_3} & \left\{ C_{p,1}^{(2)} C_j D_j \cos \left( 2p\delta_{j,1}\xi - 2\sigma_j \right) \right. \\ & + C_{p,2}^{(2)} C_3 D_3 \cos \left( 2p\delta_{j,1}\xi - 2\sigma_3 \right) \\ & + C_{p,3}^{(2)} \sqrt{C_j C_3 D_j D_3} \cos \left( 2p\delta_{j,1}\xi - \sigma_j - \sigma_3 \right) \\ & + C_{p,4}^{(2)} (C_j D_j + C_3 D_3) \\ & \left. + C_{p,5}^{(2)} \sqrt{C_j D_j C_3 D_3} \cos \left( \sigma_j - \sigma_3 \right) \right\}, \end{aligned} \quad (22) $$

<footnote n="3">³ The relevant parameter to consider is $\Gamma/G$ and the variations of $\Gamma$ are reported in G (Eqs. (37) and (38)).</footnote>
where $\delta_{j,1} = 1$ if $j = 1$ and zero otherwise. The quantities $C_{p,m}^{(n)}$ depend only on $p$ and can be obtained using the Laplace coefficients. For $p = 1$, their analytical expressions, as well as a numerical evaluation, is given in Appendix B.

The perturbation $\mathcal{H}_{1,2}$ cannot be obtained using the same procedure since the Laplace coefficients diverge for equal nominal semi-major axes and $a_{1,0} = a_{2,0} = \bar{a}$. We instead follow the method described in Robutel & Pousse (2013). We denote $\Delta = \sqrt{2 - 2 \cos \bar{\xi}}$ and find

$$ \mathcal{H}_0 = \frac{m_1}{m_0} \frac{\eta}{C_2} (\cos \xi - \Delta^{-1}), $$

$$ \mathcal{H}_{1,2}^{(2)} = \frac{m_1}{m_0} \frac{\eta}{C_2} \{A_h (C_1 D_1 + C_2 D_2) + 2 \sqrt{C_1 C_2 D_1 D_2} \mathcal{R}_c (B_h e^{i(\sigma_2 - \sigma_1)})\}, \quad (23) $$

where

$$ A_h = \frac{5 \cos 2\xi - 13 + 8 \cos \xi}{4\Delta^5} - \cos \xi \quad \text{and} \quad (24) $$

$$ B_h = e^{-2i\xi} - \frac{e^{-3i\xi} + 16e^{-2i\xi} - 26e^{-i\xi} + 9e^{i\xi}}{8\Delta^5}. $$

The final simplified Hamiltonian is then

$$ \mathcal{H} = \mathcal{H}_K + \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)}, \quad (25) $$

We denote $F_0 : \mathbb{R}^8 \to \mathbb{R}^8$ the differential system derived from Eq. (25) by the Hamilton-Jacobi equations.

## 2.2. Equilibria and linearisation in their vicinity

In this section, we study the equilibria of the resonance $p : p : p + 1$ and the dynamics in their vicinity.

### 2.2.1. Fixed points and libration centres

One of the consequences of averaging over the mean motion is that the averaged Hamiltonian (25) has equilibria, that is, points in the phase space where its gradient vanishes. The complete Hamiltonian (3) has no equilibria, however, and since the Hamiltonian (25) is supposed to model it, we are interested in the dynamics of the complete Hamiltonian at the equilibria of the model.

At a fixed point (or equilibrium) of the model, $L$ and $D_j$ are constant and so are $e_j$ and $a_j$. Similarly, the angles $\sigma_j$ and $\xi$ are constant; that is, there are constants $c_j$ such that

$$ \sigma_j = -p\lambda_2 + (p+1)\lambda_3 - \bar{\sigma}_j = c_j. \quad (26) $$

However, the secular angle $-p\lambda_2 + (p+1)\lambda_3$ and the pericentres $\bar{\sigma}_j$ are not constant at the equilibria, but they all precess with the same frequency, which we denote $\nu_3$.

The average performed in Eq. (10) is actually analogous to a first-order Lie serie expansion, and the averaging process can be seen as a periodic change of variable. Denoting $x$ and $x'$ the variables of the Hamiltonian respectively before and after the average, we have $x = e^{L_w} x'$, where $L_w = \{W, \cdot\}$ denotes the total time derivative along the trajectories of the scalar field $W$, which is constrained by the cohomological equation (Deprit 1969)

$$ \{H_K, W\} = H_P - \frac{1}{2\pi} \int_0^{2\pi} H_P \, d\xi_2. \quad (27) $$

This equation shows that, at the equilibria, $W$ is periodic of time, and so is the change of variable. This means that fixed points in the averaged model correspond to periodic trajectories in the complete system. For a quantity not invariant by rotation around the axis of the total angular momentum, however, like the secular angle $-p\lambda_2 + (p+1)\lambda_3$ or the pericentres $\bar{\sigma}_j$, a fixed point in the model corresponds in the complete Hamiltonian to a quasi-periodic motion with the two frequencies $\nu_2$ and $\nu_3$, with

$$ \nu_2 = \dot{\xi}_2 = \frac{\partial \mathcal{H}}{\partial \Gamma} \quad \text{and} \quad \nu_3 = \dot{\xi}_3 = \frac{\partial \mathcal{H}}{\partial G} = \frac{\partial \mathcal{H}}{\partial \Upsilon}, \quad (28) $$

where these quantities are evaluated at the equilibrium. More precisely, it is a periodic motion with frequency $\nu_2$ in a rotating frame following all the pericentres at frequency $\nu_3$. This result holds true for any resonance chain (e.g. Eq. (A.1) in Delisle 2017). In the rest of this work, what is known as a fixed point or equilibrium for the model is referred to as a libration centre in the complete system.

### 2.2.2. Analytical results at first order in eccentricity

Even if truncated at order 1 in eccentricity, the fixed points of the Hamiltonian (25) cannot be given analytically. Similar difficulties were met by Delisle (2017) for resonance chains with first-order resonances between non-consecutive planets. However, we show here that a further simplification allows us to obtain analytical expressions of the equilibria and of the eigenvalues of the linearised system.

To further simplify the Hamiltonian, we force a decoupling between the degree of freedom $(\xi, \Delta L)$ associated with the libration of the co-orbitals and the three other degrees of freedom $(\sigma_j, D_j)$. To this end, we first evaluate $\mathcal{H}^{(1)} = \mathcal{H}_{1,3}^{(1)} + \mathcal{H}_{2,3}^{(1)}$ at $\xi = \pi/3$. Indeed, the Hamiltonian $\mathcal{H}_K + \mathcal{H}^{(0)} + \mathcal{H}^{(1)}$ is only a perturbation of $\mathcal{H}_K + \mathcal{H}^{(0)}$ that only has equilibria at $\xi \in \{\pm\pi/3, \pi\}$, where $\xi = \pi$ is the hyperbolic (unstable) aligned configuration and $\xi = \pm\pi/3$ are the equilateral elliptic (stable) equilibria (Robutel & Pousse 2013). Both elliptic equilibria are symmetric with the same dynamics, hence we only consider $\xi = \pi/3$. Then, we replace the variable $\Delta L$ by the constant $\Delta L^*$ in the anti-diagonal term$^4$ of $\mathcal{H}_K$ in (16), where

$$ \begin{aligned} \Delta L^* &= -\frac{pC_2}{C_1 + C_2} (\Delta G + D_{1,0} + D_{2,0} + D_{3,0}) \\ &= -\frac{pC_2}{C_1 + C_2} \Delta \Upsilon^* \end{aligned} \quad (29) $$

is the value of $\Delta L$ for which $\partial \mathcal{H}_K / \partial \Delta L$ vanishes. The $D_{j,0}$ are given by Eq. (32). While the evaluation at $\xi = \pi/3$ allows analytical expressions for the position of the fixed points, the evaluation at $\Delta L = \Delta L^*$ also uncouples $(\xi, \Delta L)$ from $(\sigma_j, D_j)$ and enables analytical expressions of the eigenvalues of the linearised system in the vicinity of the fixed points. The differential system derived from $\mathcal{H}_K + \mathcal{H}^{(0)} + \mathcal{H}^{(1)}$, once these simplifications have been performed, is given in Appendix C. It vanishes when the angles are equal to

$$ \xi_0 = \frac{\pi}{3}, \quad \sigma_{1,0} = p \frac{\pi}{3} + \epsilon\pi, \quad \sigma_{2,0} = \epsilon\pi, \quad (30) $$

$$ \sigma_{3,0} = \arctan \frac{m_1 \sin p\xi_0}{m_2 + m_1 \cos p\xi_0} + (1 - \epsilon)\pi, $$

<footnote>$^4$ That is, the term proportional to $\Delta L \Delta \Upsilon$.</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Values of the eccentricities and fundamental frequencies given by the analytical results. Left: value of $e_1$ at the fixed points of the resonance chain 1:1:2 as predicted by Eq. (35). Right: values of $v/\eta$ and $v_3/\eta$ along the main branch as a function of $\delta$, for the same resonance chain, predicted by Eqs. (42) and (33). In both panels the planetary masses are $(m_1 + m_2)/2 = m_3 = 10^{-4} m_0$ and $m_2/m_1 = 10$. According to Eq. (32), $e_1 = e_2$, and for this choice of masses and resonance chain $e_3/e_1 = 0.5471$. In Table 1, which gathers the fixed points at second order in eccentricity, the three branches visible in the left plot are branches 1 (main branch), 2, and 6. The secular resonances between $v$ and $v_3$ are shown on the right, and are also easy to spot in Fig. 6.</figcaption></figure>

where

$$ \epsilon = \begin{cases} 0 & \text{if } v_3 < 0, \\ 1 & \text{if } v_3 > 0, \end{cases} \qquad (31) $$

and when the actions are equal to<fnref n="5" />

$$ \frac{C_1 C_{p,1}^{(1)2} m_1^2 p^2}{2C_3^2 m_0^2 (p+1)^2} = \left(\frac{v_3}{\eta}\right)^2 D_{1,0}, \quad C_1 D_{1,0} = C_2 D_{2,0}, \\ \frac{D_{3,0}}{D_{1,0}} = \frac{C_3 C_{p,2}^{(1)2} H^2}{C_1 C_{p,1}^{(1)2}}, \quad \Delta L_0 = 0, \qquad (32) $$

where the precession frequency of the pericentres is

$$ v_3 = -\eta K \Delta\Upsilon^*, \quad K = \frac{3p^2 C_1 C_2}{C_1 + C_2} + 3p(p+1)C_3, \qquad (33) $$

and we define the constant $H$ by

$$ H = \cos(p\xi_0 - \sigma_{3,0}) + \frac{m_2}{m_1} \cos\sigma_{3,0}. \qquad (34) $$

The unknowns of Eq. (32) are the $D_{j,0}$, and since the ratios $D_{2,0}/D_{1,0}$ and $D_{3,0}/D_{1,0}$ are known, we are reduced to the unique unknown $D_{1,0}$. The precession frequency of the pericentre, $v_3$, depends on $D_{1,0}$ (see Eqs. (29) and (33)). Denoting $C = 1 + D_{2,0}/D_{1,0} + D_{3,0}/D_{1,0}$ and performing the translation $Z = D_{1,0} + 2\Delta G/3C$, Eq. (32) is rewritten as a third-degree polynomial in $Z$:

$$ Z^3 - PZ - Q = 0, \quad \text{where} $$

$$ P = \frac{\Delta G^2}{3C^2} \quad \text{and} \quad Q = \frac{2\Delta G^3}{27C^3} + \frac{C_1 C_{p,1}^{(1)2} m_1^2 p^2}{2C_3^2 C_3^2 K^2 m_0^2 (p+1)^2}. \qquad (35) $$

The coefficients $P$ and $Q$ of this polynomial depend on the parameter $\Delta G$. There is a bifurcation between 1 and 3 real solutions when $27Q^2 - 4P^3 = 0$, that is at

$$ \Delta G = \Delta G_{\text{bif}} = - \left( \frac{27 C_1 C_{p,1}^{(1)2} p^2 C m_1^2}{8 C_3^2 K^2 m_0^2 (p+1)^2} \right)^{1/3}. \qquad (36) $$

In the rest of this work, we use the parameter

$$ \delta = \Delta G / \Delta G_{\text{bif}}, \qquad (37) $$

where

$$ \Delta G = \frac{G'}{\Gamma'} - \frac{G^*}{\Gamma^*} = \frac{G}{\Gamma} - \frac{\sum \Lambda_j^*}{\Gamma^*} = \frac{G}{\Gamma} - \sum_{j \le 3} C_j^{-1}. \qquad (38) $$

In this section, $\Gamma = \Gamma'/\Gamma^* \approx 1$ is simply evaluated at 1 and ignored, but not in Sects. 3 and 4, where tidal dissipation induces a drift in $\Gamma$, hence in $\delta$. The normalisation by $\Delta G_{\text{bif}}$ ensures that the bifurcation is at $\delta = 1$, regardless of the planetary masses.

The forced decoupling that we performed to obtain these expressions leads to results that are very similar to the second fundamental model of resonance proposed by Henrard & Lemaitre (1983). The fixed points are given by the roots of the third-degree polynomial in $Z$ (Eq. (35)), which has one or three real solutions depending on $\delta$, hence a bifurcation. The solutions of Eq. (35) are plotted in Fig. 1. For $\delta < 1$, only one elliptic equilibrium exists, called the main branch, while for $\delta \ge 1$ two other fixed points appear, one of them being hyperbolic, hence the presence of separatrices in the phase space and the formal existence of a resonance. These results come from a strong hypothesis, and we see in Sect. 2.2.3 that the topology of the Hamiltonian (25) is different (see Table 1). However, we show in Fig. 2 that these analytical expressions are accurate for low eccentricities.

<footnote n="5">⁵ $C_1 D_{1,0} = C_2 D_{2,0}$ yields $e_{1,0} = e_{2,0}$ since $e_j = \sqrt{2C_j D_j}$.</footnote>
<table><caption>**Table 1.** The 15 equilibria of the simplified Hamiltonian (25).</caption><thead><tr><th>#</th><th>100 e<sub>1</sub></th><th>100 e<sub>2</sub></th><th>100 e<sub>3</sub></th><th>σ<sub>1</sub> (°)</th><th>σ<sub>2</sub> (°)</th><th>σ<sub>3</sub> (°)</th><th>L-0.0345</th><th>ξ (°)</th><th>Nature</th><th>Domain</th></tr></thead><tbody><tr><td>1</td><td>4.449</td><td>7.490</td><td>5.878</td><td>±14.433</td><td>±22.572</td><td>±92.014</td><td>6.461e-5</td><td>±59.760</td><td>elliptic</td><td>δ ∈ ℝ</td></tr><tr><td>2</td><td>0.168</td><td>0.165</td><td>0.093</td><td>∓119.18</td><td>±179.92</td><td>±4.6606</td><td>19.34e-5</td><td>±60.003</td><td>δ-dependent</td><td>δ > 1.129</td></tr><tr><td>3</td><td>5.201</td><td>7.645</td><td>5.470</td><td>±14.216</td><td>∓15.688</td><td>±95.525</td><td>6.483e-5</td><td>±59.806</td><td>elliptic</td><td>δ > 5.997</td></tr><tr><td>4</td><td>8.344</td><td>8.646</td><td>0.305</td><td>180</td><td>0</td><td>0</td><td>6.501e-5</td><td>180</td><td>hyperbolic</td><td>δ ∈ ℝ</td></tr><tr><td>5</td><td>7.939</td><td>7.180</td><td>5.922</td><td>∓151.76</td><td>∓20.944</td><td>±92.729</td><td>6.439e-5</td><td>±179.13</td><td>hyperbolic</td><td>δ > 4.195</td></tr><tr><td>6</td><td>10.11</td><td>6.330</td><td>6.246</td><td>∓102.00</td><td>±177.92</td><td>±4.3287</td><td>7.200e-5</td><td>±62.538</td><td>hyperbolic</td><td>δ > 1.129</td></tr><tr><td>7</td><td>6.009</td><td>8.721</td><td>1.696</td><td>±6.8311</td><td>∓1.8178</td><td>±86.449</td><td>6.518e-5</td><td>±59.710</td><td>hyperbolic</td><td>δ > 5.999</td></tr><tr><td>8</td><td>0.164</td><td>0.165</td><td>0.079</td><td>0</td><td>180</td><td>0</td><td>19.34e-5</td><td>180</td><td>hyperbolic</td><td>δ > 1.082</td></tr><tr><td>9</td><td>3.708</td><td>6.584</td><td>6.887</td><td>0</td><td>180</td><td>0</td><td>7.157e-5</td><td>180</td><td>hyperbolic</td><td>δ > 1.082</td></tr></tbody></table>

**Notes.** The equilibria are found for the resonance chain 1:1:2 at $δ = 7$ using a Newton-Raphson method. Values given without decimal places are exact. The planetary masses are as in Fig. 1. Branch 2 is hyperbolic only for $5.548 ≤ δ ≤ 5.802$ and elliptic elsewhere. The entry value of $δ$ in the formal resonance (here 1.129) weakly depends on the planetary masses because of the normalisation by $ΔG_{\text{bif}}$. Branch 1 is the only elliptic branch existing for all values of $δ$ and it is the main branch introduced in Sect. 2.2.2. It corresponds to the only real solution of Eq. (35) when $δ < 1$. Branches 3, 5, and 7 do not exist at first order in eccentricity, while they exist at second order; hence, we cannot exclude that the complete Hamiltonian (3) has more libration centres, either because we did not discretise the phase space thinly enough to find them or because they do not exist at second order in eccentricity.

In the vicinity of the main branch, we linearise the differential system. We use the cartesian coordinates

$$u_j = \sqrt{2D_j} \cos \sigma_j \quad \text{and} \quad v_j = \sqrt{2D_j} \sin \sigma_j, \qquad (39)$$

and denoting $X = ^t(u_1, u_2, u_3, v_1, v_2, v_3, ΔL, ξ)$, the linearised system is

$$\begin{aligned} \frac{d\Delta X}{dt} &= \begin{pmatrix} Q_6 & 0_{6,2} \\ 0_{2,6} & Q_2 \end{pmatrix} \Delta X, \quad \text{where} \quad \Delta X = X - X_0, \\ Q_2 &= \begin{pmatrix} 0 & \frac{9m_1}{4m_0} \eta C_2^{-1} \\ -3\eta(C_1 + C_2) & 0 \end{pmatrix}, \end{aligned} \qquad (40)$$

and $X_0$ is the equilibrium value of $X$. The matrix $Q_6$ is given in Appendix D. Its characteristic polynomial reads

$$\det(\lambda \mathbb{I}_6 - Q_6) = (\lambda^2 + v_3^2)^2 (\lambda^2 + v_3^2 - 2v_3 I \sum D_{j,0}), \qquad (41)$$

where $I$ is defined in Eq. (D.1). It is interesting to note that the precession frequency of the pericentres, $±iν_3$, is an eigenvalue of the differential system (Eq. (40)). This factorisation was already noted by Pucacco (2021), who studied the resonance chain 1:2:4 of the Galilean satellites, although it was not attributed to the precession of the pericentres. The eigenvalues of $Q_2$ are $±iν$, where

$$\nu = \eta \sqrt{\frac{27 m_1 + m_2}{4 m_0}} \qquad (42)$$

is the libration frequency of the angle $ξ$ in the neighbourhood of the equilateral Lagrangian configuration (Robutel & Pousse 2013; Couturier et al. 2021). Figures 1 show that $ν_3 < 0$ for the main branch, which ensures that the roots of Eq. (41) are purely imaginary. Evaluating the eigenvalues $±iν$ and $±iν_3$ in the vicinity of the main branch shows that, at $δ ≈ -5.6$ for the planetary masses in Fig. 1, all these eigenvalues have roughly the same value, yielding a 1:1 secular resonance between the libration frequency of the co-orbitals and the precession frequency of the pericentres. Other secular resonances between $ν$ and $ν_3$ are shown in Fig. 1, and are also very visible on the stability map from Fig. 6. We show in Sect. 3.2 that the secular resonance 1:1 has important consequences for the tidal stability of the co-orbital pair.

### 2.2.3. Topology of the phase space

Limiting the work at first order in eccentricity and forcing a decoupling between $(L, ξ)$ and $(D_j, σ_j)$ gives analytical expressions of the linearised system$^6$, but at the cost of strong approximations. We develop here a Newton-Raphson-based algorithm to numerically find the equilibria of the model (25) without these approximations.

The vector field $F_0$ derived from the Hamiltonian (25) depends on the choice of the parameter $δ$ (through $ΔG$), and once a fixed point is found for a particular value of $δ$, we repeat the Newton-Raphson algorithm for slowly varying values of $δ$ in order to travel along the whole branch. We look for equilibria exploring the parallelepiped in the phase space defined by $|u_j| < 0.08$ and $|v_j| < 0.08$ ($e_j^2 = C_j(u_j^2 + v_j^2)$). We choose $L_0 = L^* = Λ_1^*/Γ^*$ as initial condition of the Newton-Raphson method for $L$ since all equilibria are expected to be close to this value, and so no discretisation is necessary along this axis. In the same way, we only choose $ξ_0 ∈ \{±π/3, π\}$.

We display in Table 1 all the equilibria that we have found for $δ = 7$, their hyperbolic or elliptic nature, and the value of $δ$ that gives birth to the branch. Due to the difficulty of exploring a thinly discretised grid in six dimensions, we may not have found all the possible equilibria. We discretised the axes $u_j$ and $v_j$ with only eight points, testing $3 × 8^6$ initial conditions, all of which converged towards 15 equilibria. Since the Hamiltonian (25) is invariant by the transformation $(ξ, σ_j) \mapsto (-ξ, -σ_j)$, fixed points with values of the angles different from 0 or $π$ have a symmetric, hence the $±$ and $∓$ signs in Table 1 (the upper sign corresponds to a fixed point and the lower sign to its symmetric). This symmetry corresponds to the invariance of the system by a rotation of angle $π$ around an axis normal to the total angular momentum.

As can be seen in Table 1, only the branches 1, 2, and 3 of fixed points can be elliptic, and thus we focus only on them in the rest of this work. In Fig. 2 we plot these branches for values of $δ$ ranging from -2 to 9. For branches 1 and 2, which are predicted by the first order in eccentricity, for comparison we also

<footnote>$^6$ Only for the equilibria and the eigenvalues; we did not obtain the eigendirections.</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Position of the elliptic branches 1, 2, and 3 (see Table 1) of equilibria of Hamiltonian (25) in the resonance 1 : 1 : 2, for $-2 \le \delta \le 9$. For each branch three curves appear, corresponding to $e_j e^{i\sigma_j}$ for $j = 1, 2, 3$. The planetary masses are the same as in Table 1. A zoomed-in image close to the origin is shown in the inset. In this area, the analytical position of these equilibria, given by Eqs. (30) and (35), is plotted by a thin grey line. The agreement is good at low values of eccentricity, but quickly worsens further from the origin. In particular, the thin grey lines are straight since $\sigma_j$ does not depend on $\delta$ in Eq. (30). Branch 1 exists for all values of $\delta$ and has all colours from yellow to dark purple, while branch 3 only exists at $\delta > 5.997$ and thus only has purple.</figcaption></figure>

plot them as given by Eqs. (30) and (35). This section shows how the analytical model is unable to locate the equilibria of the simplified Hamiltonian (25) for $e_j \gtrsim 0.005$ (see Fig. 2) and does not even give its topology for $e_j \gtrsim 0.05$ (branches 3, 5, and 7 do not exist at first order in eccentricity; see Table 1). This discrepancy between first and second order in eccentricity was already mentioned by Beaugé et al. (2006) in the case of the two-planet 1 : 2 mean motion resonance.

### 2.2.4. Comparison with the complete Hamiltonian

In this section we compare the positions of the equilibria of the secular (simplified) Hamiltonian (25) to those of the corresponding periodic orbits of the complete (full) Hamiltonian (3), which we call libration centres in Sect. 2.2.1. To this end, we develop an iterative algorithm, similar to that in Couetdic et al. (2010), to find a libration centre of the complete Hamiltonian using an equilibrium of the simplified Hamiltonian as initial condition.

We assume, close enough to a libration centre of the complete Hamiltonian (3), that the trajectories are quasi-periodic, and we write, for any complex quantity $z$ depending on these trajectories,

$$z(t) = \sum_{k \in \mathbb{Z}^6} z_k e^{ik \cdot \omega t}, \quad (43)$$

where the coordinates of $\omega = v^T(v, v_2, v_3, g_1, g_2, g_3)$ are the fundamental frequencies of the complete Hamiltonian (3). The frequencies $v$, $v_2$, and $v_3$ have approximate values given by Eqs. (42), (28), and (33), respectively.

As explained in Sect. 2.2.1, the libration centres correspond to points in the phase space where the motion is periodic$^7$ with frequency $v_2$. The description of the algorithm is as follows:

- For a choice of the parameter $\delta$, find the position of an elliptic equilibrium of the simplified Hamiltonian (25) with a Newton-Raphson method. Use it as the initial condition to integrate numerically the trajectories of the complete Hamiltonian (3);

- For a complex quantity $z$ depending on the trajectories of Eq. (3), obtain the decomposition (Eq. (43)) using a frequency analysis method (e.g. Laskar 1993);

- Identify the terms depending on frequencies other than $v_2$ and set to 0 the corresponding coefficient $z_k$.

- Proceed similarly for different quantities $z$ and evaluate them at time $t = 0$ in order to obtain a new initial condition. Restart from the first step using the new initial condition instead of the equilibrium of Eq. (25).

The process is iterated until a convergence occurs. In Fig. 3, we display the trajectories of the quantities $e_j e^{i\sigma_j}$ in the plane $(e_j \cos \sigma_j, e_j \sin \sigma_j)$ as the algorithm iterates. Isolating terms featuring frequencies other than $v_2$ is difficult, if not impossible, since the frequency analysis gives the scalars $k \cdot \omega$ but not the vector $k$, which cannot be deduced as the vector $\omega$ is unknown for the complete Hamiltonian. We can get around this difficulty because $v_2$ is much larger than the other frequencies, hence it

<footnote>$^7$ Quasi-periodic with frequencies $v_2$ and $v_3$ for a quantity not invariant by rotation around the vertical axis.</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 3. Trajectories of $e_j e^{i\sigma_j}$ in the complete Hamiltonian (3) for the six iterations needed for the algorithm to converge to the libration centre. Iteration 0 is the equilibrium (main branch) of the simplified Hamiltonian (25) at $\delta = 5$. It is rather far from the libration centre of the complete Hamiltonian (see also Fig. 4). After six iterations, the algorithm converges to the libration centre and the motion is periodic (hence the closed curves) with frequency $v_2$ (see Sect. 2.2.1). The planetary masses are the same as in Table 1 and the resonance chain is 1 : 1 : 2.</figcaption></figure>

is easy to isolate terms that depend on $v_2$ from those that do not. The implementation is thus simplified by setting to 0 the coefficients $z_k$ of the terms that do not depend on $v_2$.

This algorithm is only able to find libration centres associated with elliptic fixed points. We use it to find the branches of libration centres associated with branches 1, 2, and 3 of Table 1. We confirm the existence of a small hyperbolic zone for branch 2 in the complete Hamiltonian when the algorithm stops converging as it travels along the branch (by slowly incrementing the value of $\delta$). We plot in Fig. 4 the main branch (branch 1) of the libration centres of the complete Hamiltonian (3) and we compare it with the main branch of equilibria of the simplified Hamiltonian (25). In Fig. 5, we plot the precession frequency of the pericentres, $v_3$, for the main branch of the libration centres, and we compare it with the analytical expression (33). This figure shows that, for high values of $\delta$, $n_1/n_3$ tends towards $(p+1)/p$ (which is equal to 2 in this case). For small values of $\delta$ though, this ratio diverges from its nominal value and libration centres at small values of $\delta$ are far from the Keplerian resonance. For branch 1, we now refer to very negative values of $\delta$ as far from the resonance and to very positive values of $\delta$ as deep in the resonance. The value of $v_3$ in the complete Hamiltonian is obtained from the frequency analysis of $e^{i\xi_3}$, once the libration centre is known.

## 2.3. Stability map of the $p : p : p + 1$ resonance chain

Before taking into account tidal dissipation in the model, we study the stability of the point-mass $p : p : p + 1$ resonance

chain by constructing a dynamical map using the frequency analysis method to determine the chaoticity of a given orbit (Laskar 1990). More precisely, we study the stability of the chain along its main branch (see Fig. 1) between $\delta = 7$, deep in the resonance, and $\delta = -7$, outside the resonance.

For each value of $\delta$, we compute the position of the exact libration centre by means of the algorithm described in Sect. 2.2.4, and we choose an initial value for the angle $\xi$ between its equilibrium value (near 60°) and the value at the boundary between tadpole and horseshoe-shaped orbits (close to 24°, see Robutel & Pousse 2013). For values of $\xi_0$ close to 60°, the considered orbit is close to the main branch and it moves away for decreasing values of $\xi_0$. For all other variables, we choose as the initial condition the value at the libration centre. Every trajectory (i.e. every choice of $\delta$ and $\xi_0$) is integrated over 80 000 periods of the co-orbital planets, and for each half of the simulation the exact value of the libration frequency $v$ is extracted from the frequency analysis of $e^{i\xi}$. We obtain two values of $v$, namely $v^{(1)}$ for the first half of 40 000 periods and $v^{(2)}$ for the second half. The diffusion index, defined as (Robutel & Gabern 2006)

$$ \zeta_v = \log_{10} \left| \frac{v^{(1)} - v^{(2)}}{v^{(1)}} \right|, \qquad (44) $$

measures the degree of quasi-periodicity of the orbit. Orbits with $\zeta_v < -6$ are considered close to quasi-periodic (stable), while orbits with $\zeta_v > -2$ are very chaotic (unstable). We plot the stability map for the resonance chain 1 : 1 : 2 in Fig. 6. Secular resonances between $v$ and $v_3$, already predicted by the analytical
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Position of branch 1 of elliptic libration centres of the complete Hamiltonian (3) in the resonance 1 : 1 : 2 for $-7 \le \delta \le 7$. The planetary masses are the same as in Table 1. As a comparison, branch 1 of equilibria of the simplified Hamiltonian (25) is plotted in grey for the same range in $\delta$. For small enough values of $\delta$ the eccentricity is not too high and the agreement is good. The simplified and the complete Hamiltonian diverge when $\delta \to +\infty$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** $v_3/\eta$ as a function of $\delta$, in the complete Hamiltonian, for the resonance chain 1 : 1 : 2, along branch 1 of the libration centres in Fig. 4. The analytical expression (33) is plotted in grey for comparison. The colour gives the value of $n_1/n_3$. As expected from Eq. (33), $v_3$ is proportional to $\delta$ far from the resonance where the eccentricities are small.</figcaption></figure>

results in Fig. 1 are visible and induce chaotic motion. Overall, this stability map shows that the chain $p : p : p + 1$ is mainly stable.

### 3. Tides in the $p : p : p + 1$ resonance chain

In Sect. 2, we assumed that the bodies are point mass objects. In this section the approximation is removed and tidal dissipation due to differential and inelastic deformations of the bodies is taken into account.

#### 3.1. Extended Hamiltonian and equations of motion

Tidal contributions to the orbital evolution of the system follow a very general formulation initiated by Darwin (1880). Differential interactions between the bodies raise tidal bulges and the subsequent redistribution of mass is responsible for a

perturbation in the gravitational potential generated by body $j$ at any point $r$ in the space. This perturbation is given by (e.g. Kaula 1964)

$$V_{i,j}(\mathbf{r}) = -\kappa_2^{(j)} \frac{Gm_i}{R_j} \left(\frac{R_j}{r}\right)^3 \left(\frac{R_j}{r_i^*}\right)^3 P_2(\cos S), \quad (45)$$

where the index $i$ (resp. $j$) refers to the body responsible for the tidal bulge (resp. where the bulge is raised), $r_i$ is the position of body $i$ with respect to the barycentre of body $j$, $R_j$ is the radius of the body $j$, $\kappa_2^{(j)}$ is its second Love number, $P_2$ is the second Legendre polynomial, and $S$ is the angle between $r_i^*$ and $r$. For a body of mass $m_k$, located at $r_k$ and interacting with this bulge, the increment in potential energy is

$$U_{i,j,k}(\mathbf{r}_k) = m_k V_{i,j}(\mathbf{r}_k). \quad (46)$$

For $N = 4$ tidally interacting bodies, $N(N-1)^2 = 36$ such potentials are generated. In the case of planets orbiting a solar-type star, however, only tides raised by the star on the planets and felt by the star should be considered since they are dominant with respect to any other contribution (see Couturier et al. 2021), which means that we only consider the three contributions $U_{i,j,k}$ with $i = k = 0$ and $j \in \{1, 2, 3\}$.

The dissipation of mechanical energy inside the planets introduces a time delay $\Delta t$ between the tidal stress and the corresponding deformation. As a consequence, the tidal bulge is not aligned with the star and the subsequent torque affects the spins and orbits of the planets. For any quantity $z_j$ related to planet $j$, we write

$$z_j^*(t) = z_j(t - \Delta t_j). \quad (47)$$

For a frequency of excitation $\tau$, the quality factor $Q_j(\tau)$ (Munk & MacDonald 1960), which measures the amount of energy
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Diffusion index $\zeta_v$ as a function of $\delta$ and $\xi_0$. The planetary masses are as in Table 1 and the chain is 1 : 1 : 2. The top of the figure, at $\xi_0 \approx 60^{\circ}$, is the main branch. The blue-to-green regions are almost quasi-periodic (stable), while the red regions are chaotic (unstable). The secular resonances between the libration frequency $v$ and the precession frequency of the pericentres $v_3$, predicted by the analytical results (see Fig. 1), are visible, especially the resonance 1 : 1, which can lead to chaotic orbits for high enough values of the libration amplitude. The horseshoe-shaped orbits at the bottom are mostly chaotic, as expected, since $m_1 + m_2 = 2 \times 10^{-4}$ is close to the limit $3 \times 10^{-4}$ of their existence (Leleu et al. 2015). The main branch around the 1 : 1 resonance between $v$ and $v_3$ is tidally attractive (see Fig. 7). Systems undergoing tides entering this zone can either converge towards the top of the figure or can become completely chaotic.</figcaption></figure>

dissipated in a period $2\pi/\tau$, is related to the time delay as

$$Q_j^{-1}(\tau) = \sin(\tau\Delta t_j(\tau)) \approx \tau\Delta t_j(\tau). \qquad (48)$$

The dependency of $\Delta t_j(\tau)$ on $\tau$ is unknown, and a simple commonly used rheology consists in considering that the time delay is independent of the frequency (Mignard 1979). We adopt this tidal model in this work, reducing the rheology to the constant parameters $\kappa_2^{(j)}$ and $Q_j^{-1} = n_{j,0}\Delta t_j$.

Although tides do not preserve the total energy, the Hamiltonian formalism is extended by considering the starred variables as parameters when deriving the equations of motion. Their contribution to the Hamiltonian reads

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

where in the heliocentric reference frame

$$U_t^{(j)} = -\kappa_2^{(j)} \mathcal{G} m_0^2 \frac{R_j^5}{r_j^3 r_j^{\star 3}} P_2(\cos S), \quad T_j = \frac{\Theta_j^{\prime 2}}{2\alpha_j m_j R_j^2}, \qquad (50)$$

with

$$S = \lambda_j - \lambda_j^{\star} - (\theta_j - \theta_j^{\star}), \qquad (51)$$

and $\theta_j$ is the rotation angle of body $j$, $\omega_j^{\prime} = d\theta_j/dt$ is its rotation rate, $\Theta_j^{\prime} = \alpha_j m_j R_j^2 \omega_j^{\prime}$ is the conjugated momentum of $\theta_j$, and

$\alpha_j$ is a dimensionless structure constant depending on the state equation of body $j$ such that $\alpha_j m_j R_j^2$ is its principal moment of inertia. The transformations (9) and (11) are performed on the tidal Hamiltonian, with the normalisations $\Theta_j = \Theta_j^{\prime}/\Gamma^{\star}$ and $\mathcal{T}_j = T_j/\Gamma^{\star}$. Denoting

$$q_j = \kappa_2^{(j)} \varphi_j^5, \quad \varphi_j = \frac{R_j}{a_{j,0}}, \quad \mathcal{R}_j = C_j \Lambda_j = \frac{\bar{\Lambda}_j}{\Lambda_j^{\star}} \approx 1, \qquad (52)$$

and $\Delta\vartheta = \vartheta - \vartheta^{\star}$,

where $\vartheta$ stands for any angle, we obtain for the tidal Hamiltonian

$$\mathcal{H}_t = \sum_{j \le 3} (\mathcal{U}_t^{(j)} + \mathcal{T}_j), \qquad (53)$$

with (Couturier et al. 2021)

$$\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^{\star -6} (A_t^{(j)} + \Xi_j), \nonumber \\
\mathcal{T}_j &= \frac{C_j n_{j,0} \Theta_j^2}{2 \alpha_j \varphi_j^2}, \nonumber \\
\Xi_j &= B_t^{(j)} (\mathcal{R}_j^{-1} D_j + \mathcal{R}_j^{\star -1} D_j^{\star}) \nonumber \\
&\quad + (\mathcal{R}_j \mathcal{R}_j^{\star})^{-1/2} \mathcal{R}_e (C_t^{(j)} \sqrt{D_j D_j^{\star}} e^{i \Delta \varpi_j}), \tag{54}
\end{align}$$
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}{4} - \frac{15}{4} \cos 2(\Delta\lambda_j - \Delta\theta_j), \tag{55} \\
C_t^{(j)} &= \frac{3}{8} e^{i(\Delta\lambda_j - 2\Delta\theta_j)} + \frac{9}{4} e^{-i\Delta\lambda_j} + \frac{147}{8} e^{-i(3\Delta\lambda_j - 2\Delta\theta_j)}. \nonumber
\end{align}
$$

We note that no expansion at order 0 in the vicinity of $\Lambda_j^*$ is
performed since it loses relevant tidal dynamics (see Couturier
et al. 2021). We instead keep exact expressions in $\Lambda_j$. The dif-
ferential system is derived from the tidal Hamiltonian using the
Hamilton-Jacobi equations and considering the starred variables
as parameters. The starred variables are then expressed by a first-
order Taylor expansion in Eq. (47) (see Couturier et al. 2021 for
more details). The perturbation to the vector field, 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 \le 3} 3 \frac{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 \le 3} 3n_{j,0} \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 \\
\dot{\omega}_j &= -3n_{j,0} C_j \alpha_j^{-1} \varrho_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 their dynamics are of no interest to us, the lines $\xi_2$ and $\xi_3$ are absent from the differential system (56).

3.2. Pseudo-fixed points and linearisation in their vicinity

The total differential system that we consider for our model is
the one derived from the Hamiltonian (25), F₀, to which we now
add the tidal perturbations (Eq. (56)). We denote it F : ℝ¹³ ↪→
ℝ¹³. Here we want to find the equilibria of F and to study the
linearised dynamics in their vicinity. However, although F₀ has
equilibria (see Table 1), F has none. The five lines of Eq. (56)
corresponding to Γ̇, ĝ, ω̇₁, ω̇₂, and ω̇₃ cannot all vanish if Dⱼ ≠ 0.
Since F₀ has no equilibria at Dⱼ = 0 and does not contribute
to these five lines, we conclude that F has no equilibria. If the
planets are all synchronised, that is, if the ωⱼ all vanish<fnref n="8" />, then
Γ̇ < 0 and

$$
\dot{\Gamma} \approx - \sum_{j \le 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 \le 3} \frac{m_0}{m_j} \frac{q_j}{Q_j} e_j^2. \quad (57)
$$

We call pseudo-equilibrium of F, or pseudo-fixed
point of F, a point X ∈ ℝ¹³ such that F(X) =
⟨0, 0, 0, 0, 0, 0, 0, 0, 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 in Sect. 2.2.3.

Equation (57) shows that, on a branch of pseudo equilibria of F, the parameter Γ (and thus the parameter δ, see Eq. (38)) drifts at a speed proportional to the square of the eccentricities. This means that, with tides, the system travels along the main branch from right to left in Fig. 1 (Delisle et al. 2014) much more quickly when δ > 0 than when δ < 0 (due to high eccentricities for positive δ). 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. Since Γ is not constant at the pseudo-fixed points, but drifts at a speed given by Eq. (57), computing the eigenvalues of the linearised system makes sense only if Γ drifts slowly enough, that is, only if

$$
|\dot{\Gamma}| \ll \max_{k \le 13} |\mathcal{R}_e(\lambda_k)|, \quad (58)
$$

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

Branch 3 always has high eccentricities (see Fig. 2) and exists
only for $\delta > 5.997$. The drift in $\delta$ towards negative values is quick
at high eccentricity (see Eq. (57)), 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 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).

In Fig. 7 we plot the real parts of the eigenvalues of the linearised system associated with F, along its main branch of pseudo equilibria, which is a little perturbation of the main branch of equilibria of Eq. (25). To guarantee that the condition (58) is well respected, we limit ourselves to $\delta < 1$. This is not really a restriction since the tides ensure that this region is quickly reached. We also plot in the same figure the value of $\dot{\Gamma}$ for comparison. Only the eigenvalues that are the perturbations of Eqs. (41) and (42) are plotted. In the absence of a third planet, Couturier et al. (2021) show that the eigenvalue responsible for the exponential increase in the libration amplitude of $\xi$, and thus for the destruction of the co-orbital motion, has a real part

$$
\mathcal{R}_e(\lambda) = \frac{9}{2} \eta \frac{m_0}{m_1 + m_2} \left( \frac{m_1}{m_2} \frac{q_2}{Q_2} + \frac{m_2}{m_1} \frac{q_1}{Q_1} \right), \quad (59)
$$

and we normalise by this quantity in Fig. 7 and Appendix F.
The region $-5.64 \le \delta \le -5.47$, at the 1 : 1 secular resonance
between the libration frequency of the co-orbitals, $v$, and the fre-
quency of all the pericentres at the pseudo equilibria, $v_3$, is such
that all the eigenvalues of the system have negative real parts, and
we expect this region to be linearly stable<fnref n="9" />. We show in Sect. 4
that this is indeed the case. The linear stability is only temporary
though, since the drift in $\delta$ ensures that this region is eventually

<footnote n="8">⁸ This implies that ĝ also vanishes, since ĝ + Σ αⱼρⱼ²Cⱼ⁻¹ωⱼ = 0 by conservation of the total angular momentum.</footnote>

<footnote n="9">⁹ Even though it is slightly chaotic along the main branch for the conservative system, see the map in Fig. 6.</footnote>
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Real parts of the eigenvalues of the linearised 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 Table 1 and the resonance chain is 1 : 1 : 2. Only the four eigenvalues associated with the four degrees of freedom of the conservative system are represented. The other eigenvalues (like those associated with the rotation rates $\omega_j$) are of no interest. The real parts behave erratically at the 1 : 1 secular resonance between $\nu$ and $\nu_3$ (see Fig. 1) 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. For $\delta \le -0.104$, at least three out of four real parts are negative. In the region $\delta \le 1$ we have $|\dot{\Gamma}|/R_e(\lambda) \le 0.012$, and the criterion (58) is very well respected. Similar figures with other planetary masses are available in Appendix F.</figcaption></figure>

left. We call this the linearly stable region in the remainder of this work. The range in $\delta$ corresponding to the linearly stable region strongly depends on $m_3 / (m_1 + m_2)$. In Fig. 8 we display its position in the plane $(m_3, \delta)$.

# 4. Numerical simulations and discussions

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 also check the results drawn in Sect. 3.2 on the linearised dynamics in the vicinity of the main branch.

## 4.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 derived from Eq. (25) to which we add the tidal perturbations (Eq. (56)). The second set is a direct N-body simulation of the complete system, with the constant-$\Delta t$ model, given by the set of Eq. (E.1).

When the third planet is absent, the relevant parameters to consider to predict the destruction time of a system of two co-orbital planets are (Couturier et al. 2021)

$$ \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 \varepsilon = \frac{m_1 + m_2}{m_0}, \quad (60) $$

namely the total dissipation rate, the mass ratio, the dissipation rate ratio, and the total co-orbital 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. Due to the complexity of the dynamics of the resonance chain $p : p : p + 1$, we do

<table><caption>**Table 2.** Parameters of the two numerical simulations.</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>

**Notes.** 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 of the planets are initially synchronised 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 Sects. 2 and 3.

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, we only show the evolution of two systems for the chain 1 : 1 : 2, whose parameters are given in Table 2. Nevertheless, we performed additional simulations with different choices for the planetary masses and the initial $\delta$. The most interesting ones are presented in Appendix F, where we thoroughly discuss the influence of a larger or smaller value for $m_3$.

For the systems listed in Table 2, the chosen value of $\delta$ is such that the beginning of the simulation is at the rightmost point of the linearly stable region; these regions 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). These systems are thus expected to be initially very stable until they leave this region (due to the drift in $\delta$, see Sect. 3.2).
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** 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 linearised system associated with the vector field $F$ are computed, and the point is plotted only if all the real parts are negative. The co-orbital masses are $m_1 = m_2 = 10^{-4}$ and the tidal parameters are those of the system 0 in Table 2. 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 stabilises 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. (42) and (33), 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$ (see Appendix F for a discussion of the impact of $m_3$ on the dynamics).</figcaption></figure>

To integrate the two systems with the simplified model $F$ (Eqs. (25) and (56)), we find, for the given value of $\delta$ and the planetary masses, the position of the fixed point of the Hamiltonian (25), 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 the model with tides is very close to the fixed point of Eq. (25), and we ignore the difference. To integrate the system with the direct N-body set of Eq. (E.1), we find, for the given value of $\delta$ and the planetary masses, the position of the libration centre with the algorithm described in Sect. 2.2.4, and use it as the initial condition for the integration, again with the shift $\Delta\xi = 0.1^\circ$.

When the co-orbital planets are alone, the positivity of $R_e(\lambda)$ in Eq. (59) ensures that the system systematically reaches the horseshoe-shaped orbits and is destroyed by close encounters. In this case the time $\tau_{hs}$ needed to reach the horseshoe-shaped orbits is (Couturier et al. 2021)

$$ \tau_{hs} = \frac{\varepsilon}{9\pi\Omega} \frac{x(1+y)}{1+yx^2} \ln\left(\frac{60^\circ}{\Delta\xi}\right) T, \qquad (61) $$

where $T = 2\pi/\eta$ is the co-orbital period. Denoting $\tau_{dest}$ the co-orbital lifetime without third planet, $\tau_{hs}$ does not significantly differ from $\tau_{dest}$ for a wide range of total co-orbital mass, and

as long as $10^{-9} \le \varepsilon \le 0.005$, we have $1/2 \le \tau_{dest}/\tau_{hs} \le 2$ (Couturier et al. 2021). We can thus consider that $\tau_{hs}$ is the lifetime of the co-orbital pair, in the absence of the third planet$^{10}$. For $\Delta\xi = 0.1^\circ$, in the case of two co-orbital Earth-like planets (using the tidal parameters of Lainey 2016), we have$^{11}$

$$ \tau_{hs} = 3.771 \text{ Gyr} \left( \frac{\bar{a}}{0.04 \text{ AU}} \right)^8 \left( \frac{m_0}{m_{\odot}} \right)^{-3/2} . \qquad (62) $$

We give in Appendix G the time $\tau_{hs}$, computed from Eq. (61), for a variety of hypothetical co-orbital pairs. We normalise the time by $\tau_{hs}$ in Figs. 9, 10, 11, and in Appendix F.

## 4.2. Mechanisms of co-orbital stabilisation

In Figs. 9 and 10 the angles $\xi$ and $\sigma_j$ are plotted 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 linearised 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 dominates) guarantees a co-orbital lifetime longer than without the third planet.

In Fig. 9, 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 \text{ arcsec}$ 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 purple in Fig. 7) has been completely damped by the linearly stable region and some time is needed to pump it noticeably. Similarly, in Fig. 10, 3.5 arcsec of libration amplitude are reached at $t = 1.95 \tau_{hs}$. The early augmentation of the libration amplitude of the angles in the bottom plot is due to the fact that in the direct simulation, the linearly stable region is not exactly at the same values of $\delta$ as in the simplified model.

The linearly stable region is not the only reason why the co-orbitals in resonant chains can live longer. Another phenomenon, which we refer to as eccentricity damping stabilisation 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 stabilisation of $\xi$, due to eccentricity damping (see Fig. 11), can happen several times before horseshoe-shaped orbits are finally reached, and the system is destroyed (see Fig. 10). The explanation of the eccentricity damping stabilisation relies on the behaviour of the eccentricities. In Fig. 11, 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 \bar{\omega}, e \sin \bar{\omega})$. At time $t = t_1$ the system is still close to the fixed points, hence the quantity $e_1 e^{i\bar{\omega}t_1}$ (or any of the two other eccentricities) describes

<footnote>$^{10}$ Especially for $\varepsilon = 2 \times 10^{-4}$, for which $\tau_{dest} \approx 1.1 \tau_{hs}$.</footnote>

<footnote>$^{11}$ For the constant-Q model, the exponent is 6.5 instead of 8.</footnote>
<figure><img src="image_10.png" /><figcaption>Fig. 9. Evolution of $\xi$ and the $\sigma_j$ (in degrees) as a function of time for system 1 (see Table 2) as integrated by the simplified model F (Eqs. (25) and (56)) (top) and the direct N-body simulation (Eq. (E.1)) (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. 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 (see 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.</figcaption></figure>

a circle of small radius. Outside the linearly stable region, the eigenvalues of the linearised systems have one positive real part, and as time evolves the radius of the circle grows, while its centre (the equilibrium position of $e_1 e^{i\sigma_1}$) gets closer to the origin due to the drift in $\delta$. At time $t > t_2$ the circle $e_1 e^{i\sigma_1}$ surrounds the origin and keeps growing, as predicted by the eigenvalues, which triggers a jump in the eccentricity. On the one hand, the linearised system predicts that the circle drawn by $e_1 e^{i\sigma_1}$ grows to infinity; on the other hand, tides impose an exponential decay of the eccentricities. Indeed, the first line of Eq. (56) 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}^{-1}, \qquad (63) $$

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, which brings the system back to the vicinity of the fixed point at $t = t_4$, where the libration amplitude of the $\sigma_j$, but also of $\xi$, is small. In Figs. 9 and 10, several occurrences of the eccentricity damping stabilisation 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 stabilisation 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), 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 stabilisation 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 stabilisation of the co-orbital motion. In the absence of the third planet, Couturier et al. (2021) showed that the eccentricities are uncoupled from the co-orbital angle $\xi$. This means that the positive real part (Eq. (59)) 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. (63) still occurs, but it does not induce a stabilisation of $\xi$.

The occurrence of the eccentricity damping stabilisation 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-orbitals are destroyed before the exponential decrease in the eccentricities can save them. Deciding whether or not a given system will be saved by the eccentricity damping stabilisation requires knowing the proper modes of the linearised 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 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 stabilisation 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
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Evolution of $\xi$ and the $\sigma_j$ (in degrees) as a function of time for system 0 as integrated by the simplified model *F* (Eqs. (25) and (56)) (*top*) and by the direct *N*-body simulation (Eq. (E.1)) (*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. 6). The model reaches horseshoe-shaped orbits at $8.1 \tau_{hs}$. The grey-shaded area is the linearly stable region.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Value of $e_1$ as a function of time for system 1 in Table 2, as integrated by the simplified model *F* (Eqs. (25) and (56)) (*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 stabilisation 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. 1).</figcaption></figure>

$v$ and $v_3$ (see the stability map in Fig. 6), can lead to peculiar orbits (e.g. switching between the Lagrangian equilibria $L_4$ and $L_5$, that is, 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.

### 4.3. Discussion

It can be seen in Fig. 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 stabilisation ends up failing. This prevents the co-orbitals from living forever. We nevertheless checked that it
can work for simulations starting at $n_1/n_3 = 2.37$. Only positive values of $\delta$ allow for $n_1/n_3$ a value close to $(p+1)/p = 2$ (see Fig. 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 (see Eq. (57)) forces the system to reach the region $\delta < 0$ on a timescale much shorter than the timescale of the increase in the libration amplitude of the angles. For systems on the main branch this means that tides favour 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 & Laskar 2014).

This section shows that our simplified model (Eqs. (25) and (56)) 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 direct set of Eq. (E.1). This contrasts with the analytical work performed by Couturier et al. (2021) in the case of lone co-orbitals, where the secular model is able to quantitatively predict the outcome of the direct simulations of the complete system with less than 1% relative error (see their Table 3).

The influence of $m_3 / (m_1 + m_2)$ on the co-orbital dynamics is thoroughly discussed in Appendix F. We show that, for a large $m_3$, the linearly stable region is not efficient in stabilising the libration amplitude of $\xi$, while the eccentricity damping stabilisation is very efficient (see Fig. F.1). As $m_3$ decreases, the eccentricity damping stabilisation loses efficiency until it does not occur for very small $m_3$-values (see Fig. F.3). The linearly stable region has a maximum efficiency for $m_3 \approx 0.29(m_1 + m_2)$ (see Fig. F.2). For a small value of $m_3$, it can stabilise $\xi$ only in a tiny neighbourhood around the 1 : 1 secular resonance between $v$ and $v_3$ (see Fig. F.3). Finally, a premature destruction of the co-orbital motion (at $t < \tau_{hs}$) can occur for a very large $m_3$-value, if the system is already far from the resonance, at a $\delta$-value much lower than its value at the 1 : 1 secular resonance between $v$ and $v_3$ (see Fig. F.4).

In terms of co-orbital lifetime, the worst-case scenario occurs when the eccentricity damping stabilisation does not work and when the linearly stable region is not crossed (or is crossed while too far from the main branch). In these cases the only positive real part in the region $\delta < 0$ often has a value close to $R_e(\lambda)$ (Fig. 7), and the system reaches the horseshoe-shaped orbits at a time close to $\tau_{hs}$.

In summary, in most cases the resonance chain increases the co-orbital lifetime, but in a few cases it can also decrease it, especially when $m_3 \gg m_1 + m_2$ and $\delta$ is very negative (see Fig. F.4).

## 5. Conclusion

In this work, we have studied the dynamics of a pair of co-orbital planets in the presence of a first-order resonance with a third planet orbiting outside the co-orbitals. We have shown that for systems deep inside this resonance, many equilibria (or libration centres) exist, at least three of which are stable. The existence of a secular resonance between the libration frequency of the co-orbitals and the precession frequency of the pericentres can lead to chaotic orbits in the conservative case. However, when tides are involved, we show that this resonance stabilises the co-orbital dynamics. Another stabilisation 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, and we developed an algorithm to find them exactly. 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$, which increases the chances of a still-to-come detection of a co-orbital pair of exoplanets, since Leleu et al. (2019) have shown that co-orbital pairs are often formed within a resonance chain. While the analytical work of this paper is performed for any value of the integer $p$, figures are restricted to the chain 1 : 1 : 2 where $p = 1$. We nevertheless checked that higher values of $p$ do not impact the qualitative results.

One important contribution of this work is the discovery of a 1 : 1 secular resonance between the libration of the critical angle $\lambda_1 - \lambda_2$ and the precession of the pericentres $\boldsymbol{\omega}_j$, as well as the inherent dynamical consequences (see Figs. 6 and 7). Libration centres, quasi-periodic orbits of the unaveraged problem that generalise the equilibria of the averaged model (Sect. 2.2.1), are such that all the pericentres precess at the same frequency in their vicinity, which 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 $\boldsymbol{\omega}_j$ for the chain 1 : 2 : 3.

*Acknowledgements.* The authors thank Jean-Baptiste Delisle and Adrien Leleu for fruitful discussions and subsequent improvements of the paper. This work was supported by CFisUC (UIDB/04564/2020 and UIDP/04564/2020), GRAVITY (PTDC/FIS-AST/7002/2020), PHOBOS (POCI-01-0145-FEDER-029932), and ENGAGE SKA (POCI-01-0145-FEDER-022217), funded by COMPETE 2020 and FCT, Portugal.

## References

Beaugé, C., Michtchenko, T. A., & Ferraz-Mello, S. 2006, MNRAS, 365, 1160

Correia, A. C. M. 2009, ApJ, 704, L1

Couetdic, J., Laskar, J., Correia, A. C. M., Mayor, M., & Udry, S. 2010, A&A, 519, A10

Couturier, J., Robutel, P., & Correia, A. C. M. 2021, Celest Mech Dyn Astr, 133, 37

Cresswell, P., & Nelson, R. P. 2009, A&A, 493, 1141

Darwin, G. H. 1880, Philos. Trans. R. Soc. Lond. Ser. I, 171, 713

Delisle, J.-B. 2017, A&A, 605, A96

Delisle, J. B., & Laskar, J. 2014, A&A, 570, L7

Delisle, J.-B., Laskar, J., & Correia, A. C. M. 2014, A&A, 566, A137

Deprit, A. 1969, Celest. Mech., 1, 12

Euler, L. 1764, Novi commentarii academiae scientiarum Petropolitanae. Berlin Acad., 10, 544

Gascheau, G. 1843, C. R. Acad. Sci. Paris, 16, 393

Giuppone, C. A., Beaugé, C., Michtchenko, T. A., & Ferraz-Mello, S. 2010, MNRAS, 407, 390

Henrard, J., & Lemaitre, A. 1983, Celest. Mech., 30, 197

Kaula, W. M. 1964, Rev. Geophys. Space Phys., 2, 661

Lagrange. 1772, Œuvres complètes (Gouthier-Villars, Paris (1869))

Lainey, V. 2016, Celest. Mech. Dyn. Astron., 126, 145

Laskar, J. 1990, Icarus, 88, 266

Laskar, J. 1993, Physica D Nonlinear Phenomena, 67, 257

Laskar, J., & Robutel, P. 1995, Celest. Mech. Dyn. Astron., 62, 193

Laughlin, G., & Chambers, J. E. 2002, Astron. J., 124, 592

Leleu, A., Robutel, P., & Correia, A. C. M. 2015, A&A, 581, A128

Leleu, A., Robutel, P., & Correia, A. C. M. 2018, Celest. Mech. Dyn. Astron., 130, 24

Leleu, A., Coleman, G. A. L., & Ataiee, S. 2019, A&A, 631, A6

Mignard, F. 1979, Moon Planets, 20, 301

Munk, W. H., & MacDonald, G. J. F. 1960, The Rotation of the Earth; A Geophysical Discussion

Namouni, F. 1999, Icarus, 137, 293

Pucacco, G. 2021, Celest. Mech. Dyn. Astron., 133, 11

Robutel, P., & Gabern, F. 2006, MNRAS, 372, 1463

Robutel, P., & Pousse, A. 2013, Celest. Mech. Dyn. Astron., 117, 17
## Appendix A: Notations

For convenience, we list in Table A.1 the notations used throughout this work<fnref n="12" />.

## Appendix B: Coefficients in the expansion of the Hamiltonian

We give here the expressions of the coefficients appearing in Eqs. (21) and (22). They depend on the Laplace coefficients $b_{n/2}^m(\alpha)$ (Laskar & Robutel 1995), 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 we have

$$ 
\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) \\
&\quad + 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.1) $$

for the first order 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) \\
&\quad - 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) \\
&\quad + \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.2) $$

for the second order.

## Appendix C: Simplified differential system

We give in this appendix the expression of the differential system derived from the Hamiltonian $\mathcal{H}_K + \mathcal{H}^{(0)} + \mathcal{H}^{(1)}$, which is the Hamiltonian (25) truncated at first order in eccentricity, after the simplifications stated in Sect. 2.2.2

have been performed. We have

$$ 
\begin{aligned}
\dot{D}_1 &= -\frac{C_{p,1}^{(1)} \sqrt{2C_1 D_1} m_1 n_{3,0}}{m_0 C_3} \sin\left(p \frac{\pi}{3} - \sigma_1\right), \\
\dot{D}_2 &= -\frac{C_{p,1}^{(1)} \sqrt{2C_2 D_2} m_2 n_{3,0}}{m_0 C_3} \sin \sigma_2, \\
\dot{D}_3 &= -\frac{\sqrt{2C_3 D_3} C_{p,2}^{(1)} n_{3,0}}{m_0 C_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_0 C_1} \left(1 - \frac{1}{\Delta^3}\right), \\
\dot{\sigma}_1 &= \frac{\partial \mathcal{H}_K}{\partial \Delta \Upsilon} + \frac{C_1 C_{p,1}^{(1)} m_1 n_{3,0}}{C_3 \sqrt{2C_1 D_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_2 C_{p,1}^{(1)} m_2 n_{3,0}}{C_3 \sqrt{2C_2 D_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_3 D_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{aligned} 
\quad (C.1) $$

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

## Appendix D: Expression of the matrix $Q_6$

We give in this appendix the matrix $Q_6$ appearing in Eq. (40). 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 (D.1) $$

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 - v_3 & 0 & I s_3 s r_1 r_3 \\
0 & 0 & 0 & 0 & -v_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 - v_3 \\
-I c^2 r_1^2 + v_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 + v_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 + v_3 & -c_3 I s r_1 r_3 & 0 & -c_3 I s_3 r_3^2
\end{pmatrix} 
$$

## Appendix E: Direct complete model for tides

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{r}_j &= -\frac{\mu_j}{r_j^3} \dot{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{f_j}{\beta_j} + \sum_{k \ne j} \frac{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} 
\quad (E.1) $$

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,

<footnote n="12">¹² ὕ (qoppa) is an archaic Greek letter.</footnote>
<table><caption>Table A.1: Notations used in this paper</caption><tr><td>m<sub>0</sub>, m<sub>1</sub>, m<sub>2</sub>, m<sub>3</sub></td><td>Masses of the star, leading and trailing co-orbitals, and outermost planet</td></tr><tr><td>R<sub>1</sub>, R<sub>2</sub>, R<sub>3</sub></td><td>Radii of the leading and trailing co-orbitals and the outermost planet</td></tr><tr><td>p</td><td>Integer such that the resonance chain is p : p : p + 1</td></tr><tr><td>G, β<sub>j</sub>, μ<sub>j</sub>, μ<sub>0</sub></td><td>Gravitational constant, m<sub>0</sub>m<sub>j</sub>/(m<sub>0</sub> + m<sub>j</sub>), G(m<sub>0</sub> + m<sub>j</sub>), Gm<sub>0</sub></td></tr><tr><td>a<sub>j</sub>, e<sub>j</sub>, λ<sub>j</sub>, ω<sub>j</sub></td><td>Semi-major axis, eccentricity, mean longitude, longitude of pericentre</td></tr><tr><td>ξ, σ<sub>j</sub>, Δ, Λ<sub>j</sub></td><td>λ<sub>1</sub> - λ<sub>2</sub>, -pλ<sub>2</sub> + (p + 1)λ<sub>3</sub> - ω<sub>j</sub>, √2 - 2 cos ξ, see Eqs. (1) and (11)</td></tr><tr><td>G and Γ, C<sub>j</sub>, D<sub>j</sub></td><td>See Eqs. (9) and (11), see Eq. (15), normalised AMD (Eqs. (1), (9), and (11))</td></tr><tr><td>n<sub>j,0</sub>, a<sub>j,0</sub></td><td>Nominal mean motions, nominal semi-major axes</td></tr><tr><td>η, ¯a, ν, ν<sub>3</sub></td><td>n<sub>1,0</sub> = n<sub>2,0</sub> = (p + 1)n<sub>3,0</sub>/p, a<sub>1,0</sub> = a<sub>2,0</sub>, see Eq. (42), see Eqs. (28) and (33)</td></tr><tr><td>C<sup>(n)</sup><sub>p,m</sub>, δ, ΔG, ΔΥ*</td><td>See Appendix B, see Eq. (37), see Eq. (38), see Eq. (29)</td></tr><tr><td>κ<sub>2</sub><sup>(j)</sup>, Δt<sub>j</sub>, θ<sub>j</sub>, α<sub>j</sub></td><td>Second Love number, time-lag, rotation angle, structure coefficient</td></tr><tr><td>Q<sub>j</sub>, ϕ<sub>j</sub>, q<sub>j</sub>, R<sub>j</sub>, ω<sub>j</sub></td><td>1/(n<sub>j,0</sub>Δt<sub>j</sub>), R<sub>j</sub>/a<sub>j,0</sub>, κ<sub>2</sub><sup>(j)</sup>ϕ<sub>j</sub><sup>5</sup>, C<sub>j</sub>Λ<sub>j</sub> ≈ 1, θ̇<sub>j</sub>/n<sub>j,0</sub></td></tr><tr><td>F, F<sub>0</sub></td><td>Total differential system of the model with tides, without tides</td></tr><tr><td>R<sub>e</sub>(λ), τ<sub>hs</sub></td><td>See Eq. (59), see Eq. (61)</td></tr><tr><td>Ω, x, y, ε</td><td>See Eq. (60)</td></tr></table>

and $f_j$ is the force arising from the tidal potential energy created by the deformation of planet j (Eq. (46)):

$$f_j = -3 \frac{\kappa_2^{(j)} G m_0^2 R_j^5}{r_j^8} \mathbf{r}_j - 3 \frac{\kappa_2^{(j)} G m_0^2 R_j^5}{r_j^{10}} \Delta t_j \quad (E.2) \\ \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].$$

## Appendix F: More numerical simulations

In this appendix we present the six most interesting simulations that were not shown in Sect. 4. 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 2. We only integrate here the simplified model $F$ (Eq. (56)) (see Sect. 3.2). For each simulation, the real parts of the eigenvalues of the linearised 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$, denoted $\delta_0$, 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 normalise the real parts by $R_e(\lambda)$ (see Eq. (59)) and the times by $\tau_{hs}$ (see Eq. (61)), which is the time to reach horseshoe-shaped orbits (close to the destruction time) in the absence of a third planet. We invite the reader to have a look at Fig. 8 before reading this appendix, as $m_3$ and $\delta$ are the only varying parameters between the different simulations. 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 equilibrium.

In Fig. F.1, we have $m_3 = 8(m_1 + m_2)$, which corresponds to the value yielding the largest linearly stable region (see Fig. 8). However, for this choice of $m_3$, the linearly stable region is located at larger $\delta$-values, and so it drifts quickly (see Sect. 3.2). Furthermore, Eq. (57) shows that the drift in $\delta$ is proportional to $\sum_j m_j^{2/3}$, and so, due to the rather large $m_3$-value, the system leaves the linearly stable region after less than 0.1 $\tau_{hs}$, much more quickly than in Fig. F.2, where $m_3$ is smaller. The amplitude of $\xi$ reaches 0.07° at its lowest, not significantly smaller than its initial value of 0.1°. However, for this choice of $m_3$,

<figure><img src="image_13.png" /><figcaption>Fig. F.1: Evolution of $\xi$ and the $\sigma_j$ for $m_3 = 8(m_1 + m_2)$.</figcaption></figure>

the eccentricity damping stabilisation is very efficient, yielding a co-orbital lifetime of 8 $\tau_{hs}$.

In Fig. F.2, we have $m_3 = 0.29(m_1 + m_2)$, which is the best compromise between the width of the linearly stable region and speed of the drift in $\delta$. For this value of $m_3$, the linearly stable region splits into two distinct strips (see Fig. 8). The system stays in this region for 2 $\tau_{hs}$, 20 times longer than in Fig. F.1, and at $t = 3.2 \tau_{hs}$, $\xi$ librates with only 0.0000012° of amplitude, gaining a factor of 80 000 from the initial 0.1°. The system stays close to the main branch for 6 $\tau_{hs}$, longer than any other simulation that we performed. The eccentricity damping stabilisation is not as efficient as for $m_3 = 8(m_1 + m_2)$, but still allows the system to reach 8.5 $\tau_{hs}$ before the destruction of the co-orbital motion.

In Fig. F.3, 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$ and $\nu_3$ (the initial semi-major axes between the two simulations are very close with a tiny difference in $\delta$). Otherwise, the unique positive real part is $R_e(\lambda)$ everywhere and the destruction occurs at $\tau_{hs}$, as
<figure><img src="image_14.png" /><figcaption>Fig. F.2: Evolution of $\xi$ and the $\sigma_j$ for $m_3 = 0.29(m_1 + m_2)$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. F.3: Evolution of $\xi$ and the $\sigma_j$ for $m_3 = (m_1 + m_2)/32$.</figcaption></figure>

in the absence of the third planet, since the eccentricity damping stabilisation does not work.

In Fig. F.4, the third planet has a mass $m_3 = 19(m_1 + m_2)$ and the linearly stable region does not exist (see Fig. 8). We perform a simulation at the 1 : 1 secular resonance between $v$ and $v_3$, where it would be if it existed. In this region the positive real parts are not greater than $R_e(\lambda)$, and with the eccentricity damping stabilisation, the co-orbitals almost reach $4\tau_{hs}$. However, very far from the Keplerian resonance ($n_1/n_3 = 2.6$ at $\delta = -6$, for this choice of $m_3$), the unique positive real part is significantly greater than $R_e(\lambda)$. Furthermore, as $(L, \xi)$ and $(D_j, \sigma_j)$ are weakly coupled far from the Keplerian resonance, the eccentricity damping stabilisation is not efficient, leading to the premature destruction (here $0.8\tau_{hs}$) of the co-orbital motion. Choosing a smaller $\delta$ leads to even quicker destruction.

<figure><img src="image_16.png" /><figcaption>Fig. F.4: Evolution of $\xi$ and the $\sigma_j$ for $m_3 = 19(m_1 + m_2)$.</figcaption></figure>

## Appendix G: $\tau_{hs}$ for hypothetical co-orbital pairs

We give in Table G.1 the time $\tau_{hs}$ for hypothetical co-orbital pairs of exoplanets made up of Solar System bodies. The semi-major axis is $a = 0.04$ AU and the mass of the host star is $m_0 = m_\odot$, but $\tau_{hs}$ is easily deduced for other values using the exponents of Eq. (62). We choose $\Delta\xi = 0.1^\circ$, and again, it is straightforward to extend the results to another choice of $\Delta\xi$, since $\tau_{hs} \propto \ln(60^\circ/\Delta\xi)$ (see Eq. (61)).

<table><caption>Table G.1: $\tau_{hs}$ for some co-orbital systems.</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></tr></thead><tbody><tr><td>Earth & Earth</td><td>3.771</td><td>Earth & Moon</td><td>44.09</td></tr><tr><td>Earth & Mars</td><td>5.480</td><td>Earth & Jupiter</td><td>3.722</td></tr><tr><td>Earth & Io</td><td>2.086</td><td>Moon & Moon</td><td>50.76</td></tr><tr><td>Moon & Mars</td><td>28.16</td><td>Moon & Jupiter</td><td>50.63</td></tr><tr><td>Moon & Io</td><td>4.374</td><td>Mars & Mars</td><td>5.761</td></tr><tr><td>Mars & Jupiter</td><td>5.747</td><td>Mars & Io</td><td>2.248</td></tr><tr><td>Jupiter & Jupiter</td><td>0.7201</td><td>Jupiter & Io</td><td>2.072</td></tr><tr><td>Io & Io</td><td>2.072</td><td></td><td></td></tr></tbody></table>

In this table, the tidal parameters are those of Lainey (2016) and only the five bodies for which $\kappa_2/Q$ is well constrained are included. The close $\tau_{hs}$ between some systems is purely coincidental. It is due to the particular value of Jupiter's $\kappa_2/Q$ and to the fact that this body is much larger and much more massive than the other four.