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

HAL
open science

# Dynamics of co-orbital planets. Tides and resonance chains

Jérémy Couturier

► To cite this version:

Jérémy Couturier. Dynamics of co-orbital planets. Tides and resonance chains. Astrophysics [astro-ph]. Observatoire de Paris, 2022. English. NNT: 2022UPSLO006. tel-04197740

HAL Id: tel-04197740

https://theses.hal.science/tel-04197740

Submitted on 6 Sep 2023

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
## THÈSE DE DOCTORAT
DE L'UNIVERSITÉ PSL

Préparée à l'Institut de Mécanique Céleste et de Calcul des Éphémérides, Observatoire de Paris

# Dynamique des planètes co-orbitales
*Marées et chaînes de résonance*

# Dynamics of co-orbital planets
*Tides and resonance chains*

### Soutenue par

**Jérémy COUTURIER**

Le 21 Octobre 2022

### École doctorale n°127

**Astronomie et Astrophysique
d'Île de France**

### Spécialité

**Mécanique Céleste**

### Composition du jury :

Françoise Roques
Astronomie, Observatoire de Paris
Présidente

Alessandra Celletti
Professor, Univ. of Rome Tor Vergata
Rapportrice

Benoît Noyelles
Associate Professor, Institut UTINAM
Rapporteur

Adrien Leleu
Post-Doctorant, Université de Genève
Examinateur

Anne Lemaître
Professor, Université de Namur
Examinatrice

Philippe Robutel
Directeur de Recherche, IMCCE
Directeur de thèse

Alexandre C.M. Correia
Associate Professor, Univ. de Coimbra
Co-directeur de thèse

Gwenaël Boué
Maître de conférences, IMCCE
Invité
<figure><img src="image_3.png" /></figure>
Remerciements

Je tiens premièrement à remercier mes directeurs de thèse *Philippe Robutel* et *Alexandre C.M. Correia*. Ils ont su identifier un sujet de recherche digne de faire un beau travail de thèse, et leur encadrement optimal a permis à ces trois années de couler de source. Non seulement leurs conseils avisés m'ont permis d'élaborer une recherche de qualité, mais ils se sont aussi appliqués à rendre mes travaux visibles dans la communauté scientifique en mécanique céleste. *Philippe*, je te souhaite le meilleur en ta nouvelle qualité de directeur de l'IMCCE.

I deeply thank Françoise Roques, Alessandra Celletti, Benoît Noyelles, Adrien Leleu, Anne Lemaitre and Gwenaël Boué for being part of the jury of my oral defense. You will allow this work to officially end. Alessandra and Benoît, you deserve an additional thank for accepting to review this manuscript. Alessandra, I am grateful for the many opportunities you gave me to present my researches, be it at the CELMECH seminars, the symposium of the IAU, or even the soon-to-come CELMEC VIII at Tor Vergata. Françoise et Gwenaël, merci d'avoir géré aussi efficacement les comités de thèse.

I would like to thank Miki Nakajima and Alice Quillen for giving credit to my researches and offering me the opportunity to pursue them. I am excited to start working with you at Rochester University this fall.

Je remercie Mickaël Gastineau et Stéphane Vaillant pour leur aide sur le plan informa-
tique. Mickaël, le numérique de cette thèse doit beaucoup à tes intégrateurs "dopri8" et
"intpla", tandis que le bon déroulement de la soutenance doit à Stéphane. Merci à Jacques
Laskar d'avoir développé l'outil Trip. De manière plus générale, merci aux membres de
l'IMCCE qui permettent à ce formidable laboratoire de perdurer et d'évoluer. Merci à
Amélie Muslewski et Nevada Mendez pour l'aide dans l'organisation des missions, et à
Lusiné Amirkhanyan pour l'administratif.

Merci Antoine Petit pour ton template, la forme du manuscript te doit. Muchas
gracias a Andreas Reisenegger por la calidad de sus clases de gravitación relativista en la
PUC que me introducieron al mundo de la astrofísica. Merci à Marie-Christine Angonin
d'avoir si habilement géré mes bêtises pendant l'année de M2. Merci à Jacques Laskar de
m'avoir fait découvrir l'IMCCE en me proposant ce stage sur la théorie de la Lune, et
à Philippe Robutel, Gwenaël Boué et Hervé Manche pour leur aide dans sa complétion.
Merci à Nam Hoàng Hoài, Federico Mogavero et Léo Bernus pour les discussions entre
deux portes. Je n'oublierai pas non plus de remercier mes professeurs de maths et de
physique de MPSI et MP*, Michel Renard, Fernand Canonico, Jean Dervieux et Benoît
Gugger. La qualité de l'enseignement scientifique en classe préparatoire est trop méconnue.
Merci à Jean-Marc Lefever de m'avoir orienté vers la prépa.
Un merci considérable à Cynthia, Martin, Mathis, Mariève, Médérik et Mégane de
m'avoir accueilli si généreusement sur la glace à Iqaluit. Les résultats du Chapitre 6 vous
doivent beaucoup, surtout les bottes en Figure 6.3. Thank you so much Janine, Terry,
Dion and Lucas for such a warm welcome in your house in Tomar. I hope the construction
of the new one is going on as planned. Muchísimas gracias a todos los miembros de la
comunidad Santa Guillermina en Satipo, especialmente a Felicia, Beatriz, Jose, Daniel,
Josecito y Carmen, por los momentos maravillosos que me hicieron pasar el la selva
peruana. Gracias William por los tramites en Péru y Charo por el viaje a Lima. Gracias
Jesús por tus talentos de fotógrafo y por hacerme descubrir Lobitos. Boris, merci pour ce
voyage anticipé aux USA.

Merci à mes parents qui m'ont soutenu pendant mes études et avant, et à mes grands-
parents paternels dont la maison à été propice au travail pendant deux années de classes
préparatoires. Merci à mes grands parents maternels, à Apolline, à Baptiste, à Elly et
à Léticia pour leur gentillesse. Finalement, merci à Maxime pour ces quelques mois
formidables. J'ai hâte de voir tous les autres.
# Contents

<table><tr><td>Remerciements</td><td>i</td></tr><tr><td>1 Introduction</td><td>1</td></tr><tr><td>1.1 Definition of co-orbital motion</td><td>1</td></tr><tr><td>1.2 State of the art on co-orbital planets detection</td><td>2</td></tr><tr><td>1.3 State of the art on co-orbital planets formation</td><td>5</td></tr><tr><td>1.4 Tidal effects and motivations</td><td>5</td></tr><tr><td>1.5 Organization of the manuscript and notations</td><td>6</td></tr><tr><td>2 Main theory on celestial mechanics and dynamical systems</td><td>9</td></tr><tr><td>2.1 Framework</td><td>9</td></tr><tr><td>2.1.1 Lagrangian formalism</td><td>10</td></tr><tr><td>2.1.2 Hamiltonian formalism</td><td>11</td></tr><tr><td>2.1.3 Canonical transformation</td><td>12</td></tr><tr><td>2.1.4 Two-body problem &amp; Poincaré variables</td><td>13</td></tr><tr><td>2.2 The (n + 1)-body problem</td><td>16</td></tr><tr><td>2.2.1 Decomposition of the Hamiltonian</td><td>16</td></tr><tr><td>2.2.2 Perturbation theory : Lie serie expansion</td><td>17</td></tr><tr><td>2.3 Secular evolution of a planar non-resonant two-planet system</td><td>19</td></tr><tr><td>2.4 The differential system Ẋ = F (X)</td><td>21</td></tr><tr><td>2.4.1 Fixed points, linearization and eigenvalues</td><td>21</td></tr><tr><td>2.4.2 Eigenvalues of a perturbed matrix</td><td>22</td></tr><tr><td>2.4.3 Case where F is Hamiltonian</td><td>23</td></tr><tr><td>3 The co-orbital motion</td><td>25</td></tr><tr><td>3.1 Hamiltonian of two planar co-orbital planets</td><td>25</td></tr><tr><td>3.1.1 Expansion in the neighbourhood of the resonance</td><td>25</td></tr><tr><td>3.1.2 Averaged Hamiltonian</td><td>27</td></tr><tr><td>3.1.3 Expansion of the perturbative part</td><td>29</td></tr><tr><td>3.2 Lagrangian equilibria and linearization</td><td>30</td></tr><tr><td>3.2.1 Circular dynamics</td><td>30</td></tr><tr><td>3.2.2 Eccentric dynamics</td><td>33</td></tr><tr><td>3.3 Contribution of general relativity</td><td>35</td></tr><tr><td>3.3.1 Some intuition</td><td>35</td></tr></table>
<table>
  <tbody>
    <tr>
      <td>3.3.2 The Post-Newtonian Hamiltonian</td>
      <td>36</td>
    </tr>
    <tr>
      <td>3.3.3 Equilibria and linearization</td>
      <td>37</td>
    </tr>
    <tr>
      <td>3.4 Conclusion</td>
      <td>38</td>
    </tr>
    <tr>
      <td><b>4 Co-orbital planets in a planar first-order resonance chain</b></td>
      <td><b>39</b></td>
    </tr>
    <tr>
      <td>4.1 Hamiltonian of the p : p : p + 1 resonance chain</td>
      <td>39</td>
    </tr>
    <tr>
      <td>4.1.1 The averaged Hamiltonian</td>
      <td>39</td>
    </tr>
    <tr>
      <td>4.1.2 Expansion of the Keplerian part</td>
      <td>41</td>
    </tr>
    <tr>
      <td>4.1.3 Expansion of the perturbative part</td>
      <td>42</td>
    </tr>
    <tr>
      <td>4.2 Important theoretical results on resonance chains</td>
      <td>43</td>
    </tr>
    <tr>
      <td>4.2.1 From fixed points to libration centres</td>
      <td>43</td>
    </tr>
    <tr>
      <td>4.2.2 An algorithm to find the libration centres</td>
      <td>44</td>
    </tr>
    <tr>
      <td>4.3 Equilibria and linearization</td>
      <td>46</td>
    </tr>
    <tr>
      <td>4.3.1 Analytical results</td>
      <td>46</td>
    </tr>
    <tr>
      <td>4.3.2 Numerical results</td>
      <td>49</td>
    </tr>
    <tr>
      <td>4.4 Comparison with the unaveraged system</td>
      <td>51</td>
    </tr>
    <tr>
      <td>4.4.1 Fixed points versus libration centres</td>
      <td>51</td>
    </tr>
    <tr>
      <td>4.4.2 A stability map of the p : p : p + 1 resonance chain</td>
      <td>52</td>
    </tr>
    <tr>
      <td>4.5 Conclusion</td>
      <td>53</td>
    </tr>
    <tr>
      <td><b>5 Main theory on tidal dissipation</b></td>
      <td><b>55</b></td>
    </tr>
    <tr>
      <td>5.1 Tidal potentials</td>
      <td>55</td>
    </tr>
    <tr>
      <td>5.1.1 The perturbing potential</td>
      <td>55</td>
    </tr>
    <tr>
      <td>5.1.2 Perturbed potential and definition of the Love numbers</td>
      <td>56</td>
    </tr>
    <tr>
      <td>5.2 Simple rheologies and subsequent tidal models</td>
      <td>57</td>
    </tr>
    <tr>
      <td>5.2.1 The elastic body</td>
      <td>58</td>
    </tr>
    <tr>
      <td>5.2.2 A simple viscoelastic body : the Maxwell model</td>
      <td>60</td>
    </tr>
    <tr>
      <td>5.2.3 The constant-Δt model</td>
      <td>61</td>
    </tr>
    <tr>
      <td>5.2.4 Quality factor</td>
      <td>62</td>
    </tr>
    <tr>
      <td>5.3 The pseudo-Hamiltonian formalism</td>
      <td>63</td>
    </tr>
    <tr>
      <td><b>6 Tides in the planar three-body co-orbital problem</b></td>
      <td><b>65</b></td>
    </tr>
    <tr>
      <td>6.1 Pseudo-Hamiltonian and equations of motion</td>
      <td>65</td>
    </tr>
    <tr>
      <td>6.1.1 Expression of the pseudo-Hamiltonian</td>
      <td>65</td>
    </tr>
    <tr>
      <td>6.1.2 Equations of motions in the presence of tides</td>
      <td>67</td>
    </tr>
    <tr>
      <td>6.1.3 Conservation of the total angular momentum with tides</td>
      <td>69</td>
    </tr>
    <tr>
      <td>6.2 Equilibria and linearization</td>
      <td>69</td>
    </tr>
    <tr>
      <td>6.2.1 The eigenvalues</td>
      <td>69</td>
    </tr>
    <tr>
      <td>6.2.2 The timescales of tidal evolution</td>
      <td>72</td>
    </tr>
    <tr>
      <td>6.2.3 Looking for co-orbital exoplanets</td>
      <td>74</td>
    </tr>
    <tr>
      <td>6.3 Numerical simulations</td>
      <td>77</td>
    </tr>
    <tr>
      <td>6.3.1 Procedure</td>
      <td>77</td>
    </tr>
    <tr>
      <td>6.3.2 Comparison between the model and the complete system</td>
      <td>79</td>
    </tr>
    <tr>
      <td>6.3.3 A plateau in the energy decrease</td>
      <td>81</td>
    </tr>
    <tr>
      <td>6.4 Tides raised on the central body</td>
      <td>82</td>
    </tr>
    <tr>
      <td>6.4.1 Orders of magnitude</td>
      <td>82</td>
    </tr>
    <tr>
      <td>6.4.2 Expression of the pseudo-Hamiltonian</td>
      <td>83</td>
    </tr>
  </tbody>
</table>
<table><tbody><tr><td>6.4.3 Equilibria and eigenvalues</td><td>84</td></tr><tr><td>6.4.4 Characteristic timescales</td><td>88</td></tr><tr><td>6.5 Conclusion</td><td>89</td></tr><tr><td><b>7 Tides in the planar four-body p : p : p + 1 resonance chain</b></td><td><b>91</b></td></tr><tr><td>7.1 Pseudo-Hamiltonian and equations of motion</td><td>91</td></tr><tr><td>7.1.1 Pseudo-Hamiltonian</td><td>91</td></tr><tr><td>7.1.2 Equations of motion</td><td>93</td></tr><tr><td>7.2 Pseudo-equilibria and linearization</td><td>93</td></tr><tr><td>7.2.1 The eigenvalues</td><td>93</td></tr><tr><td>7.2.2 The linearly stable region</td><td>95</td></tr><tr><td>7.3 Numerical simulations</td><td>97</td></tr><tr><td>7.3.1 Procedure</td><td>97</td></tr><tr><td>7.3.2 Increased co-orbital life inside the resonance chain</td><td>98</td></tr><tr><td>7.3.3 The eccentricity damping stabilization mechanism</td><td>100</td></tr><tr><td>7.3.4 Influence of the mass of the third planet</td><td>103</td></tr><tr><td>7.4 Discussion and conclusion</td><td>106</td></tr><tr><td><b>8 Conclusions</b></td><td><b>109</b></td></tr><tr><td>8.1 Overview</td><td>109</td></tr><tr><td>8.2 Consequences</td><td>110</td></tr><tr><td>8.3 Perspectives</td><td>111</td></tr><tr><td><b>A Spherical harmonics and demonstration of Eq. (5.8)</b></td><td><b>113</b></td></tr><tr><td>A.1 Framework</td><td>113</td></tr><tr><td>A.2 Laplace equation in spherical coordinates</td><td>114</td></tr><tr><td>A.3 Demonstration of Eq. (5.8)</td><td>114</td></tr><tr><td><b>B Coefficients of H<sup>(4)</sup>, U<sub>t</sub><sup>(j)</sup>, H<sub>j,3</sub><sup>(1)</sup> and H<sub>j,3</sub><sup>(2)</sup></b></td><td><b>117</b></td></tr><tr><td>B.1 Coefficients of H<sup>(4)</sup></td><td>117</td></tr><tr><td>B.2 Coefficients of H<sub>j,3</sub><sup>(1)</sup> and H<sub>j,3</sub><sup>(2)</sup></td><td>117</td></tr><tr><td>B.3 Coefficients of U<sub>t</sub><sup>(j)</sup></td><td>118</td></tr><tr><td><b>C Matrices and differential systems</b></td><td><b>119</b></td></tr><tr><td>C.1 Simplified differential system</td><td>119</td></tr><tr><td>C.2 Expression of the matrix Q<sub>6</sub></td><td>120</td></tr><tr><td>C.3 Expression of the matrix Q<sub>0</sub> + Q<sub>1</sub></td><td>120</td></tr><tr><td>C.4 Complete model for tidal interactions</td><td>121</td></tr><tr><td><b>Bibliography</b></td><td><b>123</b></td></tr></tbody></table>
Introduction

1.1 Definition of co-orbital motion

In order to provide a definition of co-orbital motion, we say that two planets are co-
orbital if the difference of their mean longitude is a librating angle. More precisely, if
$\lambda_1$ and $\lambda_2$ denote the mean longitudes of two planets of a planetary system, defining
$\xi(t) = \lambda_1(t) - \lambda_2(t)$, these two planets are co-orbital if<fnref n="1" />

$$
\forall t \in \mathbb{R} \quad \xi(t) \bmod 2\pi \in ]0, 2\pi[,
\quad (1.1)
$$

where the variable *t* is the time. In the framework of the three-body problem, co-orbital
planets are in a 1 : 1 mean motion resonance. The difference ξ of their mean longitudes
generally librates around a value close to ±π/3 or equal to π. Beyond the famous collinear
and equilateral configurations described respectively by Euler (1764) and Lagrange (1772),
other types of orbits are possible in this resonance, making the co-orbital dynamics
very rich. While the Euler configurations, denoted by *L*₁, *L*₂ and *L*₃ in the restricted
three-body problem, are unstable for all mass ranges, the Lagrangian equilibria are linearly
stable, provided that the three masses satisfy the relation 27(*m*₀*m*₁ + *m*₀*m*₂ + *m*₁*m*₂) <
(*m*₀ + *m*₁ + *m*₂)² established by Gascheau (1843) for circular orbits. If we denote by *m*₀
the mass of the star, much larger than that of the planets *m*₁ and *m*₂, the Gascheau
condition is satisfied when (*m*₁ + *m*₂)/*m*₀ ≲ 1/27.

Therefore, in the planar case and for small eccentricities, when the sum of the planetary masses is smaller than about 1/27 of the stellar mass, tadpole orbits arise, allowing the difference in the mean longitude to librate around ±60°. When (m₁ + m₂)/m₀ ≲ 3 × 10⁻⁴, horseshoe-shaped orbits can arise (Laughlin and Chambers, 2002). They librate around 180° with a very large amplitude of at least 312°. For small libration amplitudes, still in the planar case, the eccentric dynamics have been studied numerically by Giuppone et al. (2010) and analytically by Robutel and Pousse (2013), showing the existence of two proper modes called Lagrange and anti-Lagrange. For moderate to large eccentricities, quasi-satellite orbits are also possible, for which the planets appear to revolve around

<footnote n="1">¹This definition excludes the quasi-satellite motion, that we do not study in this manuscript.</footnote>
each other (Giuppone et al., 2010; Pousse et al., 2017). High eccentricities give rise to topological changes in the phase space, and thus to many more exotic trajectories (Leleu et al., 2018), while the dynamics of the inclined problem is even more complex by allowing, among other things, transitions between these orbits and retrograde co-orbital motion (Namouni, 1999).

## 1.2 State of the art on co-orbital planets detection

An exoplanet is a planet orbiting a star other than the Sun. It is generally estimated that the Milky Way contains $10^{11}$ to $10^{12}$ planets, and all of them but eight are exoplanets. The era of exoplanet detection started 27 years ago, and as of August 26th, 2022, 5154 were discovered<fnref n="2" />, none co-orbital with another. Before reviewing the methods specifically used for the detection of hypothetical co-orbital exoplanets, we will focus on the detection of regular exoplanets.

Most exoplanets are located more than 10 light-years away from the Sun (only 8 out of 5154 are closer than that<fnref n="3" />), and most exoplanets have a radius less than that of Jupiter (63% of the 5154 are smaller than Jupiter<fnref n="2" />). It is likely that the proportion of exoplanets smaller than Jupiter is higher than 63%, but the detection is naturally biased towards larger exoplanets. Similarly, it is obvious that the proportion of exoplanets closer than 10 light-years is much smaller than 8/5154, since once again, the detection is strongly biased towards closer exoplanets. As a consequence, only a small minority of all exoplanets have an angular diameter larger than 0.3 milliarcseconds. Taking into account the small angular separation between the exoplanet and its host star, the feat of a direct detection is generally excluded, since the best-case scenario is equivalent to observe a fly circling a street light 14800 km away<fnref n="4" />.

As exoplanets can hardly be detected directly, their existence is generally highlighted from various kinds of perturbations they cause on their host star. If the exoplanet orbits its star in a plane almost parallel to our line of sight, then the luminosity of the star slightly decreases periodically due to eclipses caused by the planet passing in front of the star, called transits. When the radius of the planet is not too small with respect to that of the star, the S/N (signal over noise) ratio is high enough for the transit to be detected. The transit method gives the orbital period of the exoplanet with a good precision, and if the radius of the star is constrained, then that of the planet is estimated from the decrease in luminosity. Similarly, if the mass of the host star can be estimated, then the semimajor axis of the planet is deduced from Kepler's third law. However, the eccentricity and argument of pericentre generally cannot be constrained, unless the secondary eclipse, that occurs when the star passes in front of the planet, can be distinguished from the noise. The secondary eclipse is rarely detected though and only some tens of exoplanets are concerned. Some examples are given in Table A.3 of Lillo-Box et al. (2018).

Regardless of the ratio S/N, the transit method is bound to fail if the exoplanet orbits its star in a plane even slightly inclined with respect to our line of sight. In presence of an exoplanet, the star orbits the barycentre of the system star—exoplanet

<footnote n="2">²exoplanet.eu</footnote>

<footnote n="3">³Namely, Proxima Centauri b,c,d; Wolf 359 b,c & Lalande 21185 b,c.</footnote>

<footnote n="4">⁴1 cm/14800 km = 140000 km/10 ly</footnote>
periodically, and analyzing the radial velocity of the star can highlight the existence of an exoplanet. The main advantage of the radial velocity method is that it does not require the exoplanet to eclipse its host star. If the star is orbited by several exoplanets, its motion around the barycentre of the planetary system is quasiperiodic instead of periodic, and the radial velocity method is more adapted than the transit method to the detection of multi-planetary systems, as it can show the quasiperiodic motion of the star, and most planets of a planetary system might not transit. The amplitude of the motion of the star is related both to the mass *m* of the exoplanet and *i*, the inclination between its orbital plane and the plane of the sky. As a consequence, the method only allows the product *m* sin *i* to be constrained, but does not give these two quantities independently. Nevertheless, if the exoplanet detected with the radial velocity method happens to transit, then *i* = π/2 and its mass can be obtained. Contrarily to the transit method, the radial velocity method does not give constraints on the radius of the planet.

Some other techniques, still under development, can enable the detection of exoplanets, but transit and radial velocity are the two most successful. Both methods are a priori valid for the detection of co-orbital exoplanets, but difficulties inherent to the co-orbital motion arise. The transit method shows the presence of a co-orbital pair if both co-orbitals transit, but in case of moderate to large semimajor axes, even a tiny mutual inclination between the co-orbitals ensures that one of them passes above or below the star, and escapes detection. Because of the large orbital period induced by large semimajor axes, the transit method seems poorly adapted to detect pairs of co-orbital exoplanets with moderate to large semimajor axes. Small semimajor axes can allow both transits to happen despite a small mutual inclination, but as we investigate in this manuscript, close-in co-orbital pairs are subject to strong tidal dissipation that can disrupt the system.

For these reasons, the radial velocity technique seems more adapted than the transit method to locate co-orbital exoplanets. However, Leleu et al. (2017) have shown that at first order in eccentricity, a pair of co-orbital exoplanets yields the same radial velocity signature as a single exoplanet with well chosen parameters. In practice, it is hence very difficult to distinguish a pair of co-orbital exoplanets from a single exoplanet using only the radial velocity method. If both co-orbital exoplanets of masses $m_1$ and $m_2$ are at an equilateral Lagrangian equilibria, then they provoke the same perturbations on the radial velocity of the star as a single exoplanet of mass *m* located in front of $m_2$ by an angle $\varphi$ (*resp.* behind $m_1$ by an angle $\pi/3 - \varphi$), where (Leleu et al., 2017, Eq. (9))

$$m = \sqrt{m_1^2 + m_2^2 + m_1 m_2} \quad \text{and} \quad \varphi = \arctan \frac{\sqrt{3} m_1}{m_1 + 2 m_2}. \qquad (1.2)$$

As a consequence, it is impossible to differentiate the two black co-orbital exoplanets from their fictitious blue equivalent on Fig. 1.1 with the radial velocity method.

By themselves, the transit and radial velocity techniques are rather helpless on detecting co-orbital exoplanets. Combining both can nevertheless help the detection. Indeed, in the case of a single exoplanet, the transit occurs at the same time as the radial velocity reaches its mean value. In the co-orbital case though, there exists a time-shift $\Delta t$ between these two events. The time $\Delta t$ depends on the angular separation $\xi = \lambda_1 - \lambda_2$ between the co-orbitals, and denoting $\eta$ their common mean-motion, it reads (Ford and Gaudi,
<figure><img src="image_4.png" /><figcaption>Fig. 1.1 — The two black co-orbital exoplanets of masses $m_1$ and $m_2$, located at an equilateral Lagrangian equilibrium, yield on the radial velocity of the central star the same signature as the fictitious blue exoplanet. Analyzing the radial velocity of the star is thus not enough to highlight the presence of the two distinct exoplanets. If $m_1 = m_2 = m$, then the blue exoplanet of mass $\sqrt{3}m$ is at mid-distance between the co-orbitals.</figcaption></figure>

(2006; Leleu et al., 2017)

$$ \Delta t = \eta^{-1} \left( \frac{\xi}{2} + \arctan \left[ \frac{m_1 - m_2}{m_1 + m_2} \tan \frac{\xi}{2} \right] \right). \qquad (1.3) $$

If one of the co-orbital exoplanets transits, analyzing the radial velocity of the star can highlight the presence of its companion, even if it does not transit, if such a time-shift is measured. However, even a single exoplanet features a time-shift $\Delta t = -\eta^{-1}e \cos\omega$ as soon as it is eccentric with eccentricity $e$ and argument of pericentre $\omega$. Unless the quantity $e \cos\omega$ can be constrained by means of the secondary eclipse, the time-shift method is useless on detecting co-orbital exoplanets.

When the star is orbited by a single exoplanet, transits occur periodically and plotting the light curve as a function of time, with a fold at every period $T = 2\pi/\eta$, draws a vertical line, corresponding to the luminosity decrease. If the lines of pixels of the subsequent image are added, the transiting exoplanet can be detected even with a poor S/N ratio. In the case of a multi-planetary system though, the motion of each exoplanet is no longer periodic and planet-planet gravitational interactions are responsible for transit timing variations (TTV). In the presence of TTV, the line corresponding to the luminosity decrease is no longer straight but can curve according to the TTV value, giving it a river shape such as the one presented in Fig. 2 of Leleu et al. (2021). The rivers caused by the TTV induced by a co-orbital companion are very characteristic (Vokrouhlický and Nesvorný, 2014, Fig. 3) and their observation would be enough to prove the existence of a co-orbital pair. Nevertheless, in presence of large TTV, adding the lines of the river diagram does not improve the S/N ratio, and such method is not suitable to detect co-orbital exoplanets. Recently, neural networks have been shown capable of efficiently detect exoplanets in river diagrams with high TTV and low S/N (Leleu et al., 2021; Leleu et al., 2022) and could yield to the detection of a co-orbital pair of exoplanets in the future.
## 1.3 State of the art on co-orbital planets formation

The discovery of exoplanets raised the question of the existence of co-orbital planets, which are absent from the Solar System. Despite the large number of exoplanet detection, no pair or co-orbital exoplanet has been discovered so far. However, the theories of planetary formation do not prohibit the existence of co-orbital planets. Indeed, Laughlin and Chambers (2002) introduced two possible processes that can form such systems: accretion in situ or capture at a stable Lagrangian point of an already existing planet. Depending on the physical characteristics of the gas disc, the capturing process can lead to systems with a high diversity of mass ratios (Cresswell and Nelson, 2008), but also to equal-mass co-orbitals (Giuppone et al., 2012). In the in-situ scenario, different models lead to various upper limits for the mass that can form at a Lagrangian equilibrium of a giant planet. A maximum of about 0.6 Earth mass is found by Beaugé et al. (2007), while Lyra et al. (2009) obtain 5 to 15 Earth mass planets at the same location.

Once formed in the disc, the stability of the co-orbital system is not necessarily guaranteed. Beaugé et al. (2007) found that inward migration tends to slightly increase the libration amplitude of the co-orbital system, and instability during the late migrating stages with low gas friction may lead to the destruction of the system. Another study from Pierens and Raymond (2014) shows that equal-mass co-orbitals (from super-Earths to Saturns) are heavily disturbed during the gap-opening phase of their evolution. Leleu et al. (2019) studied the dynamics of a pair of migrating co-orbital planets and showed that, depending on the mass ratio, the eccentricities of the planets and the type of dissipative forces, the two planets may evolve towards the stable Lagrangian points or scatter out of the system. In the same work, Leleu et al. (2019) studied the formation of exoplanets around low mass star. They showed that when co-orbital pairs of exoplanets are formed, they are often part of a chain of mean motion resonance (see their Fig. 12).

## 1.4 Tidal effects and motivations

For systems with orbital periods less than 10 days, the planets undergo strong tidal interactions from the parent star (e.g. Correia et al., 2020), which arise from differential and inelastic deformations of the planet. In the two-body problem, the ultimate stage for tidal evolution is the synchronization of the rotational and orbital periods, alignment of the planet spin axis with the normal to the orbit, and circularization of the orbit (e.g. Hut, 1980; Adams and Bloch, 2015).

In the full N-body problem, Moeckel (2017) proved that a relative equilibrium (solid rotation of the whole configuration) is never a minimum of energy of the phase space at a given total angular momentum. Thus, applied to the three-body problem with tidal dissipation, this result implies that the Lagrangian equilibria are made unstable by tides. However, we know neither the timescale of such instability nor if the phenomenon expands to the whole phase space, and we even less know what are the consequences to the dynamics of the co-orbital configuration. Indeed, although the spin of close-in planets quickly evolves into an equilibrium configuration, the orbital evolution is much slower (e.g. Correia and Laskar, 2010), and the co-orbital configuration may survive the whole age of the system.
The absence of co-orbital exoplanet detection could be explained by the challenge that such detection represents, as we justified in Sect. 1.2. Nevertheless, in this manuscript, we rather investigate the hypothesis that tidal effects are responsible for the disruption of the co-orbital motion. Rodríguez et al. (2013) provided a numerical study of the co-orbital motion under tidal dissipation and showed that identical co-orbitals (same masses, radii and tidal parameters) are always destroyed by tides, but before Couturier et al. (2021), no analytical work had been undertaken to investigate this hypothesis. Resonance chains are at the core of current researches in celestial mechanics, and since formation models predict that co-orbital planets are often formed within one, we investigate the behaviour and stability of the co-orbital motion as part of a resonance chain.

## 1.5 Organization of the manuscript and notations

In this thesis manuscript, we provide an in-depth study of the co-orbital motion, perturbed by tides and resonance chains. The different models that we elaborate take into account effects such as the gravitational interactions between the co-orbital bodies, dissipation due to tides, secular (i.e. long-term) perturbations due to interactions with a third planet, and general relativity. The main secular effect that we are interested in is tidal dissipation. Indeed, tides do not preserve the total energy of the system and have a high potential for long-term modifications of the orbits. However, it is important that the point-mass (i.e. tideless) dynamics are well described and understood before we can properly introduce tidal dissipation. For this reason, we do not treat the case of tides before Chap. 5. In the whole manuscript, we restrict ourselves to the planar case.

The theoretical work that we perform throughout the manuscript is often analytical and sometimes semi-analytical. Numerical simulations are generally bounded to the verification of the analytical results and their domain of validity. Our work makes use of the general theory on celestial mechanics and dynamical and Hamiltonian systems. In Chap. 2, we expose the main parts of this theory, and we rederive some important results that we frequently use later in the manuscript. For convenience to the reader, we gather in Table 1.1 the notations used throughout the manuscript.

In Chap. 3, we study the dynamics of two point-mass co-orbital bodies affected by their gravitational interaction with the central body, and perturbed by their mutual interaction, that is, we study the 1:1 mean-motion resonance. Most results of Chap. 3 were first published in Robutel and Pousse (2013) and are recalled in Sect. 2 of Couturier et al. (2021). The end of Sect. 3.2.2 was first published in Couturier et al. (2021) and Sect. 3.3, where we study the perturbations due to general relativity on the co-orbital motion, is original.

Since it is not uncommon for exoplanets to be formed within a resonance chain, we expand our work in Chap. 4 by considering the case where the pair of co-orbital planets interacts with an outer third planet, in a first-order mean-motion resonance with the pair. In other words, we study the point-mass resonance chain $p: p: p+1$, where $p$ is a small integer. The results of Chap. 4 were first published in Couturier et al. (2022).

The theory of tidal dissipation could be the scope of an entire thesis manuscript, and before we introduce tides in Chaps. 6 & 7, we dedicate Chap. 5 to the main theory on tidal dissipation. We redefine there some tidal parameters that the literature often makes
use of, namely, the Love numbers and the quality factor. We also give some examples of
rheologies and their associated tidal models. In order to include tides in our model as
smoothly as possible, we explain in Sect. 5.3 how a pseudo-Hamiltonian formalism can be
developed that allows the equations of motion with tides to be derived from the Hamilton
equations.

In Chaps. 6 & 7, we include tidal dissipation to the co-orbital 1 : 1 mean-motion
resonance, and to the *p*: *p*: *p* + 1 resonance chain, respectively. We prove in Chap. 6
that the co-orbital motion is made unstable by tides, but that the third planet in the
chain *p*: *p*: *p* + 1 often makes the co-orbital pair much more stable. We give analytical
expressions of the co-orbital lifetime depending on the parameters of the system. Our
work is the first analytical proof of the instability of the Lagrangian equilibria under tides,
and the results of Chap. 6 were published in Couturier et al. (2021), while those of Chap.
7 were published in Couturier et al. (2022). Section 6.4 treats the case of tides raised on
the central body and is not yet published.
<table><caption>Table 1.1 — List of the notations used in this manuscript, ordered roughly by first apparition. A subscripted *j* means of body number *j*. Notations used right after their definition only are not included here. There exists a conflict of notation between *δ* defined by Eq. (3.14) and *δ* defined by Eq. (4.34). The conflict does not exist within a same chapter and we believe that no confusion is possible. *ξ* (varstigma), *ϑ* (qoppa) and *λ* (sampi) are archaic Greek letters.</caption><thead><tr><th>Notation</th><th>Definition</th><th>Notation</th><th>Definition</th><th>Notation</th><th>Definition</th></tr></thead><tbody><tr><td>t.</td><td>transposition</td><td>bold</td><td>vector</td><td>t</td><td>time</td></tr><tr><td>·</td><td>d/dt</td><td>○</td><td>composition</td><td>∇f</td><td>Eq. (2.2)</td></tr><tr><td>∇F</td><td>Eq. (2.4)</td><td>∇<sup>2</sup>f</td><td>Eq. (2.5)</td><td>L</td><td>Lagrangian</td></tr><tr><td>H</td><td>Hamiltonian</td><td>G</td><td>grav. constant</td><td>m<sub>j</sub></td><td>mass</td></tr><tr><td>q</td><td>coordinates</td><td>p</td><td>momenta</td><td>I</td><td>Identity matrix</td></tr><tr><td>J</td><td>Eq. (2.19)</td><td>Φ<sub>H</sub></td><td>flow of H (2.23)</td><td>β<sub>j</sub></td><td>m<sub>0</sub>m<sub>j</sub>/(m<sub>0</sub>+m<sub>j</sub>)</td></tr><tr><td>μ<sub>j</sub></td><td>G(m<sub>0</sub> + m<sub>j</sub>)</td><td>ω<sub>j</sub></td><td>long. of periapsis</td><td>e<sub>j</sub></td><td>eccentricity</td></tr><tr><td>a<sub>j</sub></td><td>semimajor axis</td><td>λ<sub>j</sub></td><td>mean longitude</td><td>Λ<sub>j</sub></td><td>(2.49) &amp; (4.10)</td></tr><tr><td>D<sub>j</sub></td><td>(2.49) &amp; (4.10)</td><td>x<sub>j</sub></td><td>Eq. (2.50)</td><td>X<sub>j</sub></td><td>√2/Λ<sub>j</sub> x<sub>j</sub></td></tr><tr><td>ι</td><td>(m<sub>1</sub> + m<sub>2</sub>)/m<sub>0</sub></td><td>L<sub>χ</sub></td><td>Eq. (2.64)</td><td>{χ,·}</td><td>L<sub>χ</sub></td></tr><tr><td>Re</td><td>real part</td><td>Jm</td><td>imaginary part</td><td>ξ</td><td>Eq. (3.7)</td></tr><tr><td>ā</td><td>Eq. (3.5)</td><td>η</td><td>Eq. (3.5)</td><td>Λ<sub>j</sub><sup>*</sup></td><td>(3.6) &amp; (4.5)</td></tr><tr><td>ξ</td><td>λ<sub>1</sub> - λ<sub>2</sub></td><td>δ</td><td>Eq. (3.14)</td><td>δ</td><td>Eq. (4.34)</td></tr><tr><td>J</td><td>(3.18) &amp; (3.20)</td><td>J<sub>2</sub></td><td>(3.18) &amp; (3.20)</td><td>R<sub>j</sub></td><td>Eq. (3.20)</td></tr><tr><td>f<sub>j</sub></td><td>Eq. (3.21)</td><td>m</td><td>√m<sub>1</sub>m<sub>2</sub></td><td>Δ</td><td>√2 - 2 cos ξ</td></tr><tr><td>A<sub>h</sub></td><td>Eq. (3.30)</td><td>B<sub>h</sub></td><td>Eq. (3.30)</td><td>ν</td><td>Eq. (3.32)</td></tr><tr><td>g<sub>1</sub></td><td>Eq. (3.42)</td><td>g<sub>2</sub></td><td>Eq. (3.42)</td><td>c</td><td>speed of light</td></tr><tr><td>n<sub>j</sub><sup>(0)</sup></td><td>Eq. (4.4)</td><td>ξ<sub>2,3</sub></td><td>Eq. (4.7)</td><td>σ<sub>j</sub></td><td>Eq. (4.7)</td></tr><tr><td>L,G,Γ</td><td>below Eq. (4.7)</td><td>Γ<sup>*</sup></td><td>Eq. (4.9)</td><td>C<sub>j</sub></td><td>Eq. (4.14)</td></tr><tr><td>p</td><td>small integer</td><td>C<sub>p,m</sub><sup>(n)</sup></td><td>Append. (B.2)</td><td>F<sub>0</sub></td><td>below Eq. (4.21)</td></tr><tr><td>ν<sub>2</sub></td><td>Eq. (4.24)</td><td>ν<sub>3</sub></td><td>Eq. (4.30)</td><td>ΔG</td><td>Eq. (4.35)</td></tr><tr><td>B</td><td>perturbed body</td><td>W</td><td>perturbing pot.</td><td>V</td><td>perturbed pot.</td></tr><tr><td>ζ</td><td>Eq. (5.16)</td><td>k<sub>lm</sub></td><td>Love number</td><td>h</td><td>Eq. (5.16)</td></tr><tr><td>μ</td><td>shear modulus</td><td>ρ</td><td>density</td><td>g</td><td>surface gravity</td></tr><tr><td>κ<sub>2</sub></td><td>Eq. (5.36)</td><td>Q</td><td>Sect. 5.2.4</td><td>Δt</td><td>time lag</td></tr><tr><td>.*</td><td>Eq. (6.5)</td><td>R<sub>j</sub></td><td>radius</td><td>Δζ</td><td>ζ - ζ<sup>*</sup></td></tr><tr><td>ϑ<sub>j</sub></td><td>(6.7) &amp; (7.2)</td><td>q<sub>j</sub></td><td>κ<sub>2</sub><sup>(j)</sup> q<sub>j</sub><sup>5</sup></td><td>Ω</td><td>q<sub>1</sub>/Q<sub>1</sub> + q<sub>2</sub>/Q<sub>2</sub></td></tr><tr><td>C<sub>j</sub></td><td>mom. of inertia</td><td>α<sub>j</sub></td><td>C<sub>j</sub>/(m<sub>j</sub>R<sub>j</sub><sup>2</sup>)</td><td>ϑ<sub>j</sub></td><td>Eq. (6.14)</td></tr><tr><td>A<sub>t</sub><sup>(j)</sup></td><td>Eq. (6.10)</td><td>B<sub>t</sub><sup>(j)</sup></td><td>Eq. (6.10)</td><td>C<sub>t</sub><sup>(j)</sup></td><td>Eq. (6.10)</td></tr><tr><td>λ<sub>j</sub></td><td>Eq. (6.28)</td><td>λ</td><td>Eq. (6.28)</td><td>λ<sub>L</sub></td><td>Eq. (6.28)</td></tr><tr><td>λ<sub>AL</sub></td><td>Eq. (6.28)</td><td>τ<sub>L</sub></td><td>Eq. (6.29)</td><td>τ<sub>AL</sub></td><td>Eq. (6.29)</td></tr><tr><td>x</td><td>m<sub>1</sub>/m<sub>2</sub></td><td>y</td><td>q<sub>2</sub>Q<sub>1</sub>/(q<sub>1</sub>Q<sub>2</sub>)</td><td>τ<sub>lib</sub></td><td>Eq. (6.35)</td></tr><tr><td>Δξ</td><td>above Eq. (6.37)</td><td>τ<sub>hs</sub></td><td>Eq. (6.40)</td><td>M<sub>⊕</sub></td><td>Earth mass</td></tr><tr><td>M<sub>☉</sub></td><td>Solar mass</td><td>R<sub>⊕</sub></td><td>Earth radius</td><td>κ<sub>p</sub><sup>q</sup></td><td>Eq. (6.59)</td></tr><tr><td>λ</td><td>Eq. (6.61)</td><td>ā<sub>c</sub></td><td>Eq. (6.62)</td><td>ā<sub>s</sub></td><td>Eq. (6.63)</td></tr><tr><td>F</td><td>Sect. 7.1.2</td><td>ω<sub>j</sub></td><td>below Eq. (7.8)</td><td>δ<sub>0</sub></td><td>δ(t = 0)</td></tr></tbody></table>
# Main theory on celestial mechanics and dynamical systems

This chapter establishes the main mathematical framework on celestial mechanics and dynamical systems, that we use throughout the rest of this manuscript.

## 2.1 Framework

We consider for this chapter a system of $n + 1$ point masses $m_j$ ($0 \le j \le n$), moving under the effect of their mutual gravitational interactions. Their positions and speeds in an inertial reference frame are given by $\mathbf{r}_j$ and $\dot{\mathbf{r}}_j$, respectively. In this manuscript, $^t$ denotes the transpose operator and an upper dot denotes $d/dt$. If $p$ is an integer and $f$ a differentiable scalar field

$$f : \mathbb{R}^p \to \mathbb{R} \\ ^t(x_1, \dots, x_p) \mapsto f(^t(x_1, \dots, x_p)) \tag{2.1}$$

then $\nabla f : \mathbb{R}^p \to \mathbb{R}^p$ is defined as

$$\nabla f = {}^t \left( \frac{\partial f}{\partial x_1}, \dots, \frac{\partial f}{\partial x_p} \right) \tag{2.2}$$

Writing $\boldsymbol{x}$ for $^t(x_1, \dots, x_p)$, we will also use the notation $\nabla f = \partial f / \partial \boldsymbol{x}$. If $\boldsymbol{x}_k$ represents the subvector $^t(x_1, \dots, x_k)$, where $k \le p$, then $\nabla_{\boldsymbol{x}_k} f : \mathbb{R}^p \to \mathbb{R}^k$ is defined as

$$\nabla_{\boldsymbol{x}_k} f = {}^t \left( \frac{\partial f}{\partial x_1}, \dots, \frac{\partial f}{\partial x_k} \right) \tag{2.3}$$

and once again, we will use the notation $\nabla_{\boldsymbol{x}_k} f = \partial f / \partial \boldsymbol{x}_k$. If $F : \mathbb{R}^p \to \mathbb{R}^p$ denotes a differentiable vector field, then, writing $F = (F_1, \dots, F_p)$, where the $F_i$ are scalar fields
like in Eq. (2.1), the notation $\nabla F$ is the Jacobian matrix of $F$

$$ \nabla F = \begin{pmatrix} \frac{\partial F_1}{\partial x_1} & \cdots & \frac{\partial F_1}{\partial x_p} \\ \vdots & \ddots & \vdots \\ \frac{\partial F_p}{\partial x_1} & \cdots & \frac{\partial F_p}{\partial x_p} \end{pmatrix}. \qquad (2.4) $$

Finally, if $f$ is twice differentiable, we denote $\nabla^2 f$ the Hessian matrix of $f$, that is

$$ \nabla^2 f = \begin{pmatrix} \frac{\partial^2 f}{\partial x_1 \partial x_1} & \cdots & \frac{\partial^2 f}{\partial x_1 \partial x_p} \\ \vdots & \ddots & \vdots \\ \frac{\partial^2 f}{\partial x_p \partial x_1} & \cdots & \frac{\partial^2 f}{\partial x_p \partial x_p} \end{pmatrix}. \qquad (2.5) $$

### 2.1.1 Lagrangian formalism

We denote $\boldsymbol{q} = {}^t(\boldsymbol{r}_0, \dots, \boldsymbol{r}_n)$ and $\dot{\boldsymbol{q}} = {}^t(\dot{\boldsymbol{r}}_0, \dots, \dot{\boldsymbol{r}}_n)$. The Lagrangian of the system is defined as $\mathcal{L} = \mathcal{T} - \mathcal{U}$, where $\mathcal{T}$ and $\mathcal{U}$ are the kinetic and potential energies, respectively. They read

$$ \mathcal{T} = \sum_{j=0}^{n} \frac{1}{2} m_j \dot{\boldsymbol{r}}_j \cdot \boldsymbol{r}_j \quad \text{and} \quad \mathcal{U} = - \sum_{i=0}^{n} \sum_{j=i+1}^{n} \frac{\mathcal{G} m_i m_j}{|\boldsymbol{r}_i - \boldsymbol{r}_j|}, \qquad (2.6) $$

where $\mathcal{G}$ is the gravitational constant. We define $\boldsymbol{x} = {}^t(\boldsymbol{q}, \dot{\boldsymbol{q}})$. The state $\boldsymbol{x}(t)$ of the system at time $t$ is given by a point in $\mathbb{R}^{6n+6}$, and its evolution from $t = t_0$ to $t = t_1$ is a path

$$ \gamma: [t_0, t_1] \to \mathbb{R}^{6n+6} \\ t \mapsto \boldsymbol{x}(t). \qquad (2.7) $$

The action $S(\gamma)$ of the path $\gamma$ is defined as

$$ S(\gamma) = \int_{t_0}^{t_1} \mathcal{L}(\boldsymbol{q}, \dot{\boldsymbol{q}}) dt, \qquad (2.8) $$

and, according to the least action principle, the path actually followed by the system extremizes $S$, that is

$$ \frac{\delta S}{\delta \gamma}(\boldsymbol{h}) := \lim_{\epsilon \to 0} \frac{S(\gamma + \epsilon \boldsymbol{h}) - S(\gamma)}{\epsilon} = 0, \qquad (2.9) $$

for every loop $\boldsymbol{h}: [t_0, t_1] \to \mathbb{R}^{6n+6}$ such that $\boldsymbol{h}(t_0) = \boldsymbol{h}(t_1) = 0$. Writing $\boldsymbol{h} = (h, \dot{h}) \in ([t_0, t_1] \to \mathbb{R}^{3n+3})^2$, we obtain

$$ \frac{\delta S}{\delta \gamma}(\boldsymbol{h}) = \int_{t_0}^{t_1} \frac{\partial \mathcal{L}}{\partial \boldsymbol{q}} \cdot h(t) + \frac{\partial \mathcal{L}}{\partial \dot{\boldsymbol{q}}} \cdot \dot{h}(t) dt, \qquad (2.10) $$

and taking advantage of $h(t_0) = h(t_1) = 0$, an integration by part yields

$$ \frac{\delta S}{\delta \gamma}(\boldsymbol{h}) = \int_{t_0}^{t_1} \left( \frac{\partial \mathcal{L}}{\partial \boldsymbol{q}} - \frac{d}{dt} \frac{\partial \mathcal{L}}{\partial \dot{\boldsymbol{q}}} \right) \cdot h(t) dt = 0. \qquad (2.11) $$
Since Eq. (2.11) is verified for any choice of the loop $h$, the least action principle is equivalent to the Euler-Lagrange equation

$$ \frac{\partial \mathcal{L}}{\partial \boldsymbol{q}} - \frac{d}{dt} \frac{\partial \mathcal{L}}{\partial \dot{\boldsymbol{q}}} = 0. \qquad (2.12) $$

Substituting Eq. (2.6) into Eq. (2.12) yields Newton's law of gravitation

$$ \ddot{\boldsymbol{r}}_k = - \sum_{j \neq k} \frac{\mathcal{G} m_j}{|\boldsymbol{r}_k - \boldsymbol{r}_j|^3} (\boldsymbol{r}_k - \boldsymbol{r}_j), \qquad (2.13) $$

but the Euler-Lagrange equation applies to many domains. As an example, the equation of geodesics in general relativity can be retrieved from Eq. (2.12) using<fnref n="1" /> $\mathcal{L} = {}^t\mathcal{u}gu$, where $g$ is the spacetime metric.

## 2.1.2 Hamiltonian formalism

Given a Lagrangian $\mathcal{L}(\boldsymbol{q}, \dot{\boldsymbol{q}})$, the generalized momenta $p_j$ are defined as

$$ \boldsymbol{p} = {}^t(\boldsymbol{p}_0, \dots, \boldsymbol{p}_n) = \frac{\partial \mathcal{L}}{\partial \dot{\boldsymbol{q}}}, \qquad (2.14) $$

while the Hamiltonian $\mathcal{H}$ is given by the Legendre transformation

$$ \mathcal{H} = \boldsymbol{p} \cdot \dot{\boldsymbol{q}} - \mathcal{L}. \qquad (2.15) $$

Substituting Eq. (2.14) into the differential of Eq. (2.15) yields

$$ d\mathcal{H} = \dot{\boldsymbol{q}} \cdot d\boldsymbol{p} - \frac{\partial \mathcal{L}}{\partial \boldsymbol{q}} \cdot d\boldsymbol{q}, \qquad (2.16) $$

showing that $\mathcal{H}(\boldsymbol{p}, \boldsymbol{q})$ depends only on $\boldsymbol{p}$ and $\boldsymbol{q}$, and giving

$$ \dot{\boldsymbol{q}} = \frac{\partial \mathcal{H}}{\partial \boldsymbol{p}} \quad \text{and} \quad \frac{\partial \mathcal{H}}{\partial \boldsymbol{q}} = -\frac{\partial \mathcal{L}}{\partial \boldsymbol{q}}. \qquad (2.17) $$

Combining together Eqs. (2.12), (2.14) & (2.17), we obtain the Hamilton equations

$$ \dot{\boldsymbol{q}} = \frac{\partial \mathcal{H}}{\partial \boldsymbol{p}} \quad \text{and} \quad \dot{\boldsymbol{p}} = -\frac{\partial \mathcal{H}}{\partial \boldsymbol{q}}. \qquad (2.18) $$

The vector $X = {}^t(\boldsymbol{p}, \boldsymbol{q})$ belongs to a 6 ($n+1$) dimensional space, called phase space of the system. The number of degrees of freedom of the system, that we denote $N$, is defined as half the dimension of the phase space. The Hamilton equations are often written under the compact form

$$ \dot{X} = \mathcal{J} \nabla \mathcal{H}(X) = \begin{pmatrix} 0 & -\mathbb{I}_N \\ \mathbb{I}_N & 0 \end{pmatrix} \begin{pmatrix} \partial \mathcal{H} / \partial \boldsymbol{p} \\ \partial \mathcal{H} / \partial \boldsymbol{q} \end{pmatrix}. \qquad (2.19) $$

<footnote n="1">¹with $\boldsymbol{u} = {}^t(\frac{dt}{d\tau}, \frac{dx}{d\tau}, \frac{dy}{d\tau}, \frac{dz}{d\tau})$, $\tau$ the proper time and $d/d\tau$ instead of $d/dt$ in Eq. (2.12).</footnote>
### 2.1.3 Canonical transformation

When studying the dynamics of a planetary system, arbitrary changes of variables do not preserve, a priori, the form of the Hamilton equations. A transformation

$$ \Psi : \mathbb{R}^{2N} \to \mathbb{R}^{2N} \\ X \mapsto Y \tag{2.20} $$

is said to be canonical if it verifies

$$ \dot{X} = \mathcal{J} \nabla \mathcal{H}(X) \quad \Rightarrow \quad \dot{Y} = \mathcal{J} \nabla \tilde{\mathcal{H}}(Y), \tag{2.21} $$

where $\tilde{\mathcal{H}} = \mathcal{H} \circ \Psi^{-1}$. That is, $\Psi$ is canonical if it preserves the form of Eqs. (2.18). We denote $\mathcal{M} = \nabla \Psi = (\partial Y_i / \partial X_j)_{1 \le i,j \le 2N}$ the Jacobian matrix of $\Psi$. Verifying that $\dot{Y} = \mathcal{M} \dot{X}$ and $\nabla \mathcal{H}(X) = {}^t \mathcal{M} \nabla \tilde{\mathcal{H}}(Y)$, and combining with Eq. (2.19), we obtain

$$ \dot{Y} = \mathcal{M} \mathcal{J}^t \mathcal{M} \nabla \tilde{\mathcal{H}}(Y). \tag{2.22} $$

A matrix $\mathcal{Q}$ is said to be symplectic if $\mathcal{Q} \mathcal{J}^t \mathcal{Q} = \mathcal{J}$. As a consequence, a transformation is canonical if, and only if, its Jacobian is symplectic.

The flow of $\mathcal{H}(X)$ is defined as the application

$$ \Phi_{\mathcal{H}} : \mathbb{R} \times \mathbb{R}^{2N} \to \mathbb{R}^{2N} \\ (t, X_0) \mapsto X(t), \tag{2.23} $$

where $X(t)$ is the unique solution to the Cauchy problem $\dot{X} = \mathcal{J} \nabla \mathcal{H}(X)$ with initial condition $X(0) = X_0$. We can give a sufficient condition of canonicity using the notion of flow. A transformation $\Psi : \mathbb{R}^{2N} \to \mathbb{R}^{2N}$ is canonical if there exists a scalar $t \in \mathbb{R}$ and a Hamiltonian $\chi : \mathbb{R}^{2N} \to \mathbb{R}$, called generator of $\Psi$, such that (Morbidelli, 2002, Sect. 1.6)

$$ \Psi = \Phi_{\chi}(t, \cdot). \tag{2.24} $$

This result can be proven in two steps :

* We first prove that the Jacobian matrix $\mathcal{M}(t) = \nabla \Phi_{\chi}(t, X)$ of the transformation verifies the differential equation $\dot{\mathcal{M}} = \mathcal{J} \mathcal{S} \mathcal{M}$, where $\mathcal{S}$ is a symmetrical matrix.

* We then prove that $\mathcal{M}$ is symplectic, which ends the demonstration, since we already proved that a transformation is canonical if, and only if, its Jacobian is symplectic.

#### First step

The new set of variables $Y(t) = \Psi(X(t)) = \Phi_{\chi}(t, X)$ is given by the flow of $\chi$, thus we have

$$ \frac{d}{dt} \Phi_{\chi}(t, X) = \mathcal{J} \nabla \chi(\Phi_{\chi}(t, X)). \tag{2.25} $$

If we evaluate Eq. (2.25) at $X + dX$, where $dX$ is a small displacement of the initial condition, a first-order Taylor expansion gives

$$ \frac{d}{dt} \left( \Phi_{\chi}(t, X) + \mathcal{M}(t) dX \right) = \mathcal{J} \nabla \chi \left( \Phi_{\chi}(t, X) + \mathcal{M}(t) dX \right). \tag{2.26} $$

<footnote n="2">²Or equivalently, if ${}^t \mathcal{Q} \mathcal{J} \mathcal{Q} = \mathcal{J}$, since $\mathcal{J}^2 = -\mathbb{I}_{2N}$.</footnote>
Another Taylor expansion in the right-hand side yields $\dot{M}(t)dX = \mathbb{J}S(t)M(t)dX$, where

$$S(t) = \nabla^2 \chi(\Phi_{\chi}(t, X)) \qquad (2.27)$$

is the Hessian of $\chi$ (Eq. (2.5)), evaluated at $\Phi_{\chi}(t, X)$, which is a symmetrical matrix. This last equation being valid for any value of $dX$, we obtain

$$\dot{M}(t) = \mathbb{J}S(t)M(t). \qquad (2.28)$$

## Second step

The condition stating that $M$ is symplectic is $^tM\mathbb{J}M = \mathbb{J}$. We have

$$\frac{d}{dt} (^tM\mathbb{J}M) = ^t\dot{M}\mathbb{J}M + ^tM\dot{\mathbb{J}}M = ^tM (^tS^t\mathbb{J}\mathbb{J} + \mathbb{J}^2S)M = 0, \qquad (2.29)$$

where we used Eq. (2.28) for the second equality and $^t\mathbb{J}\mathbb{J} = \mathbb{I} = -\mathbb{J}^2$ for the third one. Integrating both sides of Eq. (2.29) with respect to time and using $M(0) = \mathbb{I}$, we obtain

$$^tM\mathbb{J}M = \mathbb{J}, \qquad (2.30)$$

and $M(t)$ is symplectic, and $X \mapsto \Phi_{\chi}(t, X)$ is canonical.

This sufficient condition for canonicity will be very useful in Sect. 2.2.2, where we introduce Lie series. We can give another sufficient condition for canonicity, very useful in practice. If $\Psi: \mathbb{R}^{2N} \to \mathbb{R}^{2N}$ is linear and takes the form

$$\Psi: \begin{pmatrix} \mathbf{p} \\ \mathbf{q} \end{pmatrix} \mapsto \left( \begin{array}{c|c} \mathcal{A} & 0 \\ 0 & {}^{t}\mathcal{A}^{-1} \end{array} \right) \begin{pmatrix} \mathbf{p} \\ \mathbf{q} \end{pmatrix} = \nabla \Psi \begin{pmatrix} \mathbf{p} \\ \mathbf{q} \end{pmatrix}, \qquad (2.31)$$

where $\mathcal{A}$ is any $N \times N$ non-degenerate matrix, then it is canonical. This is readily seen by checking that the Jacobian $\nabla\Psi$ is symplectic. We end this section by showing that a Hamiltonian $\mathcal{H}(\mathbf{p}, \mathbf{q})$ is conserved along its own flow. Indeed, we have

$$\frac{d\mathcal{H}}{dt} = \frac{\partial \mathcal{H}}{\partial \mathbf{p}} \cdot \dot{\mathbf{p}} + \frac{\partial \mathcal{H}}{\partial \mathbf{q}} \cdot \dot{\mathbf{q}} = \dot{\mathbf{q}} \cdot \dot{\mathbf{p}} - \dot{\mathbf{p}} \cdot \dot{\mathbf{q}} = 0. \qquad (2.32)$$

## 2.1.4 Two-body problem & Poincaré variables

### The two-body problem

Combining Eqs. (2.6), (2.14) & (2.15), the Hamiltonian of the two-body problem can be written

$$\mathcal{H}_{2BP}(\mathbf{p}, \mathbf{q}) = \frac{\mathbf{p}_0 \cdot \mathbf{p}_0}{2m_0} + \frac{\mathbf{p}_1 \cdot \mathbf{p}_1}{2m_1} - \frac{\mathcal{G}m_0m_1}{|\mathbf{q}_0 - \mathbf{q}_1|}. \qquad (2.33)$$

This Hamiltonian has 6 degrees of freedom, and 3 of them can be eliminated by performing the canonical$^3$ transformation

$$(\mathbf{q}_0, \mathbf{q}_1; \mathbf{p}_0, \mathbf{p}_1) \mapsto (\mathbf{r}_0, \mathbf{r}; \tilde{\mathbf{r}}_0, \tilde{\mathbf{r}}) = (\mathbf{q}_0, \mathbf{q}_1 - \mathbf{q}_0; \mathbf{p}_0 + \mathbf{p}_1, \mathbf{p}_1). \qquad (2.34)$$

<footnote>$^3$It has the form (2.31).</footnote>
The new Hamiltonian does not depend on $\mathbf{r}_0$, and given Eqs. (2.18), this means that $\tilde{\mathbf{r}}_0 = \mathbf{p}_0 + \mathbf{p}_1$ is constant. Terms depending only on $\tilde{\mathbf{r}}_0$ can therefore be removed from the Hamiltonian, since their gradient is zero. We obtain

$$ \mathcal{H}_{2BP}(\tilde{\mathbf{r}}, \mathbf{r}) = \frac{\tilde{\mathbf{r}}^2}{2\beta} - \frac{\tilde{\mathbf{r}}_0 \cdot \tilde{\mathbf{r}}}{m_0} - \frac{Gm_0m_1}{r}, \qquad (2.35) $$

where $\beta = m_0m_1/(m_0 + m_1)$. The Hamiltonian can be further simplified by the translation$^4$ $\hat{\mathbf{r}} = \tilde{\mathbf{r}} - \beta\tilde{\mathbf{r}}_0/m_0$. This yields

$$ \mathcal{H}_{2BP}(\hat{\mathbf{r}}, \mathbf{r}) = \frac{\hat{\mathbf{r}}^2}{2\beta} - \frac{Gm_0m_1}{r}. \qquad (2.36) $$

The Hamiltonian $\mathcal{H}_{2BP}$ in Eq. (2.36) still has 3 degrees of freedom, but we can additionally eliminate two of them. We first notice that the angular momentum $\mathbf{G} = \mathbf{r} \times \hat{\mathbf{r}}$ is conserved

$$ \frac{d\mathbf{G}}{dt} = \dot{\mathbf{r}} \times \hat{\mathbf{r}} + \mathbf{r} \times \frac{d\hat{\mathbf{r}}}{dt} = \frac{1}{\beta}\hat{\mathbf{r}} \times \hat{\mathbf{r}} - \frac{Gm_0m_1}{r^3}\mathbf{r} \times \mathbf{r} = 0. \qquad (2.37) $$

By definition of the cross product, $\mathbf{G}$ is perpendicular to both $\mathbf{r}$ and $\hat{\mathbf{r}}$. As a consequence, the motion occurs in a plane perpendicular to $\mathbf{G}$ containing the origin, and the Hamiltonian $\mathcal{H}_{2BP}$ is thus reduced to two degrees of freedom. We write

$$ \mathbf{r} = {}^t(u, v) \quad \text{and} \quad \hat{\mathbf{r}} = {}^t(\hat{u}, \hat{v}), \qquad (2.38) $$

where $u$ and $v$ (resp. $\hat{u}$ and $\hat{v}$) are the cartesian coordinates of $\mathbf{r}$ (resp. $\hat{\mathbf{r}}$) in the orbital plane. The transformation to the polar coordinates $(u, v; \hat{u}, \hat{v}) \mapsto (r, \varphi; \tilde{r}, \Phi)$, where

$$ u = r \cos \varphi, \quad v = r \sin \varphi; \quad \tilde{r} = \frac{u\hat{u} + v\hat{v}}{r} \quad \text{and} \quad \Phi = u\hat{v} - v\hat{u}, \qquad (2.39) $$

is canonical (Meyer and Hall, 2009, Sect. 7.4). By definition of the cross product, $\Phi$ is the norm of the total angular momentum

$$ \Phi = |\mathbf{r} \times \hat{\mathbf{r}}| = |\mathbf{G}| = G. \qquad (2.40) $$

Since $G = \Phi$ is conserved, the Hamiltonian $\mathcal{H}_{2BP}$ (Eq. (2.36)) in the new coordinates does not depend on $\varphi$ and reads

$$ \mathcal{H}_{2BP}(\tilde{r}, r) = \frac{1}{2\beta} \left( \tilde{r}^2 + \frac{G^2}{r^2} \right) - \frac{Gm_0m_1}{r}. \qquad (2.41) $$

We are reduced to one degree of freedom. The canonical coordinates $(r, \varphi; \tilde{r}, G)$, where $\varphi$ is the true anomaly, are the Hill coordinates (Hill, 1913; Laskar, 2017).

<footnote>$^4$The Jacobian is identity, and the identity is symplectic, so it is a canonical transformation.</footnote>
Poincaré variables

The form of Hamilton Eqs. (2.18) shows that the flow of a Hamiltonian $\mathcal{H} = \mathcal{H}(\mathbf{p})$ depending only on the momenta is trivial. It is given by

$$ \mathbf{p}(t) = \mathbf{p}(0) \quad \text{and} \quad \mathbf{q}(t) = \mathbf{q}(0) + t\omega_0 = \mathbf{q}(0) + t\frac{\partial \mathcal{H}}{\partial \mathbf{p}}. \qquad (2.42) $$

We now transform the Hamiltonian of the two-body problem so that it only depends on the momentum. We know since Newton that the orbit is a conic<fnref n="5" />

$$ r(\varphi) = \frac{a(1-e)(1+e)}{1+e\cos\varphi}, \qquad (2.43) $$

where $a$ is the semimajor axis of the conic and $e$ its eccentricity. At the periapsis ($\varphi = 0$) and apoapsis ($\varphi = \pi$), the value of $r$ is

$$ r_p = a(1-e) \quad \text{and} \quad r_a = a(1+e). \qquad (2.44) $$

Furthermore, at these points in the trajectory, $r$ reaches an extremum and $\tilde{r} = \beta\dot{r}$ is zero. The conservation of the Hamiltonian along the trajectory (Eq. (2.32)) allows us to write $\mathcal{H}_{2BP}(0, r_p) = \mathcal{H}_{2BP}(0, r_a)$, which yields

$$ G = \beta\sqrt{\mu a(1-e^2)}, \quad \text{where} \quad \mu = \mathcal{G}(m_0 + m_1). \qquad (2.45) $$

Evaluating the Hamiltonian at the periapsis, or at the apoapsis, gives

$$ \mathcal{H}_{2BP}(a) = -\frac{\mu\beta}{2a}, \qquad (2.46) $$

and $\mathcal{H}_{2BP}$ depends only on $a$. However, the transformation to the elliptic elements $(\hat{\mathbf{r}}; \mathbf{r}) \mapsto (a, e, i; M, \omega, \Omega)$, where $M$ is the mean anomaly<fnref n="6" />, $i$ the inclination of the orbit, $\omega$ the argument of the pericentre and $\Omega$ the longitude of the ascending node, is not canonical. Delaunay (1861) & Andoyer (1923) solve the problem by introducing the canonical transformation $(\hat{\mathbf{r}}; \mathbf{r}) \mapsto (L, G, H; M, \omega, \Omega)$, where Delaunay variables read

$$ \begin{aligned} \Lambda &= \beta \sqrt{\mu a} & M, \\ G &= \beta \sqrt{\mu a (1 - e^2)} & \omega, \\ H &= G \cos i & \Omega. \end{aligned} \qquad (2.47) $$

See Laskar (2017) for a proof of canonicity. In Delaunay coordinates, the Hamiltonian takes the expression

$$ \mathcal{H}_{2BP}(\Lambda) = -\frac{\mu^2 \beta^3}{2\Lambda^2}. \qquad (2.48) $$

This Hamiltonian has one degree of freedom and depends only on the momenta. All the Delaunay coordinates are thus constant, except $M(t) = nt = t\partial\mathcal{H}_{2BP}/\partial\Lambda = \mu^2\beta^3t/\Lambda^3$ which evolves linearly with time. To avoid the equality of momenta when either $e$ or $i$

<footnote n="5">⁵This comes from the conservation of the eccentricity vector $e = (m_0 + m_1)\hat{\mathbf{r}} \times \mathbf{G}/(\mathcal{G}m_0^2m_1^2) - \mathbf{r}/r$.</footnote>

<footnote n="6">⁶Defined as $M(t) = 2\pi t/T$, where $T$ is the period of the orbit.</footnote>
is zero, we introduce Poincaré polar coordinates through a canonical linear change of variable of the form (2.31)

$$ \begin{array}{l} \Lambda = \beta \sqrt{\mu a} \\ D = \Lambda - G = \Lambda (1 - \sqrt{1 - e^2}) \\ U = G - H = G (1 - \cos i) \end{array} \quad \left| \begin{array}{l} \lambda = M + \varpi, \\ -\varpi = -\omega - \Omega, \\ -\Omega, \end{array} \right. \qquad (2.49) $$

where $\lambda$ and $\varpi$ are the mean longitude and longitude of the pericentre, respectively. In the rest of this manuscript, we only consider the planar case, and we do not use the variables $U$ and $\Omega$. We can also adopt complex Poincaré variables with the canonical transformation

$$ (D; -\varpi) \mapsto (x = \sqrt{D} \exp(i\varpi); -i\bar{x}), \qquad (2.50) $$

where the upper bar denotes the complex conjugated. In Chap. 3, we start from Poincaré complex variables, while we rather use the polar Poincaré variables in Chap. 4.

## 2.2 The $(n+1)$-body problem

### 2.2.1 Decomposition of the Hamiltonian

We combine Eqs (2.6), (2.14) & (2.15) to give the Hamiltonian of the $(n+1)$-body problem under the form

$$ \mathcal{H}(\tilde{\mathbf{u}}, \mathbf{u}) = \sum_{j=0}^{n} \frac{\tilde{\mathbf{u}}_j \cdot \tilde{\mathbf{u}}_j}{2m_j} - \sum_{j=0}^{n} \sum_{i=j+1}^{n} \frac{\mathcal{G}m_i m_j}{|\mathbf{u}_i - \mathbf{u}_j|}. \qquad (2.51) $$

As for the two-body problem, three degrees of freedom can be eliminated by identifying body $j$ with respect to body 0, for $j \neq 0$. We thus perform the transformation

$$ \mathbf{u}_j \mapsto \mathbf{r}_j = \mathbf{u}_j - \mathbf{u}_0 \quad \text{if } j \neq 0, \quad \text{and} \quad \mathbf{u}_0 \mapsto \mathbf{r}_0 = \mathbf{u}_0. \qquad (2.52) $$

The transformation is made canonical by transforming the momenta $\tilde{\mathbf{u}}_j \mapsto \tilde{\mathbf{r}}_j$ according to Eq. (2.31). The resulting Hamiltonian does not depend on $\mathbf{r}_0$ and $\tilde{\mathbf{r}}_0$ is constant. It contains terms of the form $\tilde{\mathbf{r}}_0 \cdot \tilde{\mathbf{r}}_j$, and those can be eliminated by an extension of the translation of momenta that we performed for the two body problem. We obtain the $3n$ degrees of freedom Hamiltonian (Laskar and Robutel, 1995)

$$ \mathcal{H}(\tilde{\mathbf{r}}, \mathbf{r}) = \sum_{j=1}^{n} \left( \frac{\tilde{\mathbf{r}}_j^2}{2\beta_j} - \frac{\mathcal{G}m_0 m_j}{r_j} \right) + \sum_{j=1}^{n} \sum_{i=j+1}^{n} \left( \frac{\tilde{\mathbf{r}}_i \cdot \tilde{\mathbf{r}}_j}{m_0} - \frac{\mathcal{G}m_i m_j}{|\mathbf{r}_i - \mathbf{r}_j|} \right), \qquad (2.53) $$

where $\beta_j = m_0 m_j / (m_0 + m_j)$. The first term of the Hamiltonian is a sum of two-body Hamiltonians like (2.36), and according to Eq. (2.48), Eq. (2.53) can be written in Poincaré polar coordinates as

$$ \mathcal{H}(\Lambda, D; \lambda, -\varpi) = \mathcal{H}_K(\Lambda) + i\mathcal{H}_P(\Lambda, D; \lambda, -\varpi), \qquad (2.54) $$

where we restricted ourselves to the planar case$^7$,

$$ \mathcal{H}_K(\Lambda) = - \sum_{j=1}^{n} \frac{\beta_j^3 \mu_j^2}{2\Lambda_j^2}, \qquad (2.55) $$

<footnote>$^7$And so, the Hamiltonian (2.54) has only $2n$ degrees of freedom.</footnote>
and $\mu_j = G(m_0 + m_j)$. In most practical problems of planetary dynamics, one body is much more massive than the others, and $m_j \ll m_0$, for $j \neq 0$. The Hamiltonian (2.54) is thus a sum of Kepler problems, perturbed by $i\mathcal{H}_P$. The Keplerian Hamiltonian $\mathcal{H}_K$ contains, for instance, star-planet interactions while $i\mathcal{H}_P$ contains the planet-planet interactions. The dimensionless parameter $i$ is defined as

$$i = \frac{m_1 + m_2 + \dots + m_n}{m_0}, \qquad (2.56)$$

and is introduced to emphasize the fact that the perturbation is much smaller than the sum of Kepler problems. Contrary to the Keplerian part, the perturbation $i\mathcal{H}_P$ cannot be given analytically as a function of the Poincaré coordinates, but Laskar and Robutel (1995) give the expansion

$$i\mathcal{H}_P = \sum_{i=1}^{n} \sum_{j=i+1}^{n} i\mathcal{H}_{i,j}, \qquad (2.57)$$

where

$$i\mathcal{H}_{i,j} = \sum_{\mathbf{k} \in \mathbb{Z}^2} \left( \sum_{\mathbf{q} \in \mathbb{N}^4} \Xi(\Lambda_i, \Lambda_j) X_i^{q_1} X_j^{q_2} \bar{X}_i^{q_3} \bar{X}_j^{q_4} \right) e^{i(k_1 \lambda_i + k_2 \lambda_j)} \qquad (2.58)$$

corresponds to the interaction between planet $i$ and planet $j$. In this expression, $X_j = \sqrt{2/\Lambda_j} x_j = e_j e^{i\varpi_j} + \mathcal{O}(e_j^3)$, where $x_j$ is the complex Poincaré coordinate introduced by Eq. (2.50). Most of the coefficients $\Xi(\Lambda_i, \Lambda_j)$ are zero. Indeed, the conservation of the total angular momentum implies that the Hamiltonian is invariant by the transformation$^8$ $(\lambda_i, \lambda_j, \varpi_i, \varpi_j) \mapsto (\lambda_i + \vartheta, \lambda_j + \vartheta, \varpi_i + \vartheta, \varpi_j + \vartheta)$, where $\vartheta$ is any constant angle. Injecting this transformation into Eq. (2.58) shows that, for a non-zero $\Xi$, the tuples $\mathbf{k} = (k_1, k_2) \in \mathbb{Z}^2$ and $\mathbf{q} = (q_1, q_2, q_3, q_4) \in \mathbb{N}^4$ comply with the d'Alembert rule

$$k_1 + k_2 + q_1 + q_2 - q_3 - q_4 = 0. \qquad (2.59)$$

### 2.2.2 Perturbation theory : Lie series expansion

We suppose in this subsection that we are facing a $N$ degrees of freedom Hamiltonian similar to that of Eq. (2.54)

$$\mathcal{H}(\mathbf{I}, \boldsymbol{\theta}) = \mathcal{H}_0(\mathbf{I}) + i\mathcal{H}_1(\mathbf{I}, \boldsymbol{\theta}), \qquad (2.60)$$

where $\mathbf{I} = (I_1, \dots, I_N)$ is the vector of momenta (also called actions), while $\boldsymbol{\theta} = (\theta_1, \dots, \theta_N)$ are the associated angular coordinates. Since $i$ is much smaller than unity, this Hamiltonian depends almost only on the momenta, and if $i\mathcal{H}_1$ were to be neglected, it would be trivially integrated by Eq. (2.42). We assume that among the angles of $\boldsymbol{\theta}$, one of them, say $\theta_1$, is circulating at high frequency and we would like to lose one degree of freedom and simplify the dynamics by averaging$^9$ the Hamiltonian over it. In other words, we want to find a canonical transformation

$$\Psi : \mathbb{R}^{2N} \to \mathbb{R}^{2N} \\ (\mathbf{I}, \boldsymbol{\theta}) \mapsto (\mathbf{J}, \boldsymbol{\varphi}), \qquad (2.61)$$

<footnote>$^8$This is an occurrence of Noether's theorem.</footnote>

<footnote>$^9$The procedure is readily generalized when several (or all) angles are averaged.</footnote>
close to the identity, such that

$$ \tilde{\mathcal{H}}(\mathbf{J}, \varphi) = \mathcal{H}(\mathbf{I}, \theta) = \tilde{\mathcal{H}}_0(\mathbf{J}) + \iota \tilde{\mathcal{H}}_1(\mathbf{J}, \varphi_2, \dots, \varphi_N) + \iota^2 \tilde{\mathcal{H}}_2(\mathbf{J}, \varphi), \quad (2.62) $$

where $\tilde{\mathcal{H}} = \mathcal{H} \circ \Psi^{-1}$. The dependency on $\varphi_1$ is thus pushed at the second order in $\iota$, and at first order, the subsequent Hamiltonian has only $N-1$ degrees of freedom left. To ensure a canonical transformation, we choose (Eq. (2.24))

$$ \Psi = \Phi_{\iota\chi}(-1, \cdot), \quad \text{or equivalently} \quad \Psi^{-1} = \Phi_{\iota\chi}(1, \cdot), \quad (2.63) $$

where $\iota\chi$ is a generator that we will constrain to achieve the form (2.62). Once again, $\iota$ is written to emphasize the order of magnitude of the generator. Such an order of magnitude makes sure that $\Psi$ is $\iota$-close to the identity. We define the Lie derivative of $\chi(\mathbf{J}, \varphi)$ as the time derivative along the flow generated by $\chi$

$$ L_{\chi} := \left. \frac{\partial \chi}{\partial \mathbf{J}} \cdot \frac{\partial}{\partial \varphi} - \frac{\partial \chi}{\partial \varphi} \cdot \frac{\partial}{\partial \mathbf{J}} \right|_{\chi} = \dot{\varphi} \cdot \frac{\partial}{\partial \varphi} + \dot{\mathbf{J}} \cdot \frac{\partial}{\partial \mathbf{J}} = \frac{d}{dt} \Big|_{\chi}. \quad (2.64) $$

The Lie derivative often appears written $L_{\chi} = \{\chi, \cdot\}$ in the literature$^{10}$. We can now give for $\tilde{\mathcal{H}}$ the expression

$$ \tilde{\mathcal{H}} = \mathcal{H} \circ \Psi^{-1} = \mathcal{H} \circ \Phi_{\iota\chi}(1, \cdot) = \sum_{j=0}^{+\infty} \frac{L_{\iota\chi}^{(j)}}{j!} (\mathcal{H}) = e^{\iota L_{\chi}}(\mathcal{H}), \quad (2.65) $$

where $L_{\iota\chi}^{(j)}$ denotes the $j$th iteration of the Lie derivative, or the $j$th time-derivative along the flow of $\iota\chi$. Let us expand Eq. (2.65) at second order in $\iota$. We obtain

$$ \tilde{\mathcal{H}}_0 + \iota \tilde{\mathcal{H}}_1 + \iota^2 \tilde{\mathcal{H}}_2 = \mathcal{H}_0 + \iota (\mathcal{H}_1 + \{\chi, \mathcal{H}_0\}) + \iota^2 (\{\chi, \mathcal{H}_1\} + \frac{1}{2} \{\chi, \{\chi, \mathcal{H}_0\}\}), \quad (2.66) $$

where, given Eq. (2.65), the different Poisson brackets are evaluated in the new variables. This yields the following equation, said cohomological, to constrain the generator

$$ \mathcal{H}_1 - \tilde{\mathcal{H}}_1 = \{\mathcal{H}_0, \chi\}, \quad (2.67) $$

where we impose on $\tilde{\mathcal{H}}_1$ to be the average of $\mathcal{H}_1$ over $\varphi_1$

$$ \tilde{\mathcal{H}}_1 = \frac{1}{2\pi} \int_0^{2\pi} \mathcal{H}_1(\mathbf{J}, \varphi) d\varphi_1. \quad (2.68) $$

The cohomological Eq. (2.67) allows $\chi$ to be constrained, and $\tilde{\mathcal{H}}_2$ to be subsequently determined with Eq. (2.66). The dependency on $\varphi_1$ can be pushed to even larger orders of $\iota$ by repeating the procedure. The Hamiltonian $\tilde{\mathcal{H}}_0 + \iota \tilde{\mathcal{H}}_1$ becomes the new principal part (formerly $\mathcal{H}_0$), while $\iota^2 \tilde{\mathcal{H}}_2$ becomes the new perturbative part (formerly $\iota \mathcal{H}_1$), and a new iteration of the method is performed to push the dependency on $\varphi_1$ to the third order on $\iota$. The process could, a priori, be pushed up to infinity in order to find a canonical transformation that completely removes the dependency on $\varphi_1$, but in general, the sequence of remainders $\iota^n \tilde{\mathcal{H}}_n$ diverges beyond a certain order and there exists an optimal order on $\iota$.

<footnote>$^{10}$Where the operator $\{\cdot, \cdot\}$ is called Poisson bracket.</footnote>
Often enough in practice, one will only be interested in pushing the dependency on $\varphi_1$ up to $\iota^2$ and to consider $\iota^2 \tilde{\mathcal{H}}_2$ as a negligible remainder. In that case, the knowledge of $\chi$ is only useful to retrieve the old variables as a function of the new ones, a transformation seldomly performed since the new variables are more desirable than the old ones. Indeed, the canonical transformation induced by the method is quasiperiodic and close to identity, and the new variables give the tendency of the trajectory, where the fast oscillations due to $\varphi_1$ have been removed. This means that in most cases, a first-order Lie series expansion reduces to a mere average of the Hamiltonian over the fast-circulating angles.

When needed, the cohomological Eq. (2.67) is solved in the Fourier domain. We write

$$ \chi = \sum_{k \in \mathbb{Z}^N \setminus \{0\}} \chi_k e^{ik \cdot \varphi}, \quad \text{and} \quad \mathcal{H}_1 - \tilde{\mathcal{H}}_1 = \sum_{k \in \mathbb{Z}^N \setminus \{0\}} h_k e^{ik \cdot \varphi}, \qquad (2.69) $$

and denoting $\omega_0 = \partial \mathcal{H}_0 / \partial \mathbf{J}$ the frequency vector, the solutions of Eq. (2.67) are

$$ \chi_k = \frac{h_k}{ik \cdot \omega_0}. \qquad (2.70) $$

For a non-zero $h_k$, it is important that $\omega_0$ respects the Diophantine-like condition $k \cdot \omega_0 > \iota$, so that the order of magnitude of $\iota\chi$ is $\iota$ and the transformation is $\iota$-close to the identity. If, for instance, the system features a mean-motion resonance $\xi_{i,j} = p_i \lambda_i + p_j \lambda_j \approx 0$, where $p_i$ and $p_j$ are small integers, then a linear canonical transformation has first to be applied to the Hamiltonian in order for $\xi_{i,j}$ to appear explicitly. Then, the Lie expansion is performed without averaging upon $\xi_{i,j}$. We use the results of this section at several occasions throughout the manuscript, and especially in Sect. 4.2.1 where we explain the notion of *libration centre*.

## 2.3 Secular evolution of a planar non-resonant two-planet system

The remaining chapters of this work will focus on the dynamics of planetary systems where mean motion resonances between planets exist. The 1:1 resonance will be the core of Chaps. 3 & 6, while the chain of resonance 1:1:2, and more generally $p:p:p+1$, will be studied in Chaps. 4 & 7. In this section, we focus on the simple case of a planar two-planet system, far from any mean motion resonance, limiting ourselves to the second order in eccentricity. This is known as the Laplace-Lagrange approximation. The star and two planets have masses $m_0, m_1$ and $m_2$, respectively. According to Eq. (2.54), the Hamiltonian of the problem in Poincaré canonical complex variables (Eq. (2.50)) reads

$$ \mathcal{H} = - \sum_{j \in \{1,2\}} \frac{\mu_j^2 \beta_j^3}{2\Lambda_j^2} + i \mathcal{H}_P(X_1, X_2, \bar{X}_1, \bar{X}_2, \lambda, \Lambda), \qquad (2.71) $$

where $X_j = \sqrt{2/\Lambda_j} x_j$. We are mostly interested in the secular (i.e. long-term) evolution of the system, and so, we can eliminate two out of the four degrees of freedom by averaging over the mean longitudes $\lambda_1$ and $\lambda_2$. As we proved in Sect. 2.2.2, this amounts in performing a canonical transformation on the Hamiltonian that sends the dependency on
these variables in a remainder of order $\iota^2$, that we neglect. For convenience, and as we do not plan to perform the inverse transformation$^{11}$, we do not change the name of the variables. In Eq. (2.58), averaging over $\lambda_1$ and $\lambda_2$ is equivalent to discarding every terms such that $k \neq 0$. Taking into account the d'Alembert rule (2.59), we can deduce that the expansion of the perturbative part, truncated at second order in eccentricity, only features terms of the form $X_j \bar{X}_j$ or $X_j \bar{X}_k$, for $j, k \in \{1, 2\}$, and we obtain the decomposition

$$ \mathcal{H} = - \sum_{j \in \{1,2\}} \frac{\mu_j^2 \beta_j^3}{2\Lambda_j^2} + \frac{m_1 n_2 \Lambda_2}{m_0} \left( B_0 + B_1 (X_1 \bar{X}_1 + X_2 \bar{X}_2) + B_2 (X_1 \bar{X}_2 + X_2 \bar{X}_1) \right), \quad (2.72) $$

where $n_2$ is the mean motion of the second planet, defined just below Eq. (2.48). The dimensionless coefficients $B_0, B_1$ and $B_2$ depend on the $\Lambda_j$, that are now constant, since the dependency on the $\lambda_j$ has been lost. As a consequence, the Keplerian part and the factor $B_0$ are now constant, their gradient is zero and they can be removed from the Hamiltonian without changing the equations of motion. We thus have the final form of the two-degrees of freedom Laplace-Lagrange Hamiltonian

$$ \mathcal{H}_{\text{LL}}(X_1, X_2, \bar{X}_1, \bar{X}_2) = \frac{m_1 n_2 \Lambda_2}{m_0} \left( B_1 (X_1 \bar{X}_1 + X_2 \bar{X}_2) + B_2 (X_1 \bar{X}_2 + X_2 \bar{X}_1) \right). \quad (2.73) $$

Laskar and Robutel (1995) (top of page 216) give the coefficients $B_1$ and $B_2$ as a function of the Laplace coefficients

$$ \begin{aligned} B_1 &= -\frac{1}{8}\alpha b_{3/2}^{(1)}(\alpha), \\ B_2 &= \frac{1}{4}(1+\alpha^2)b_{3/2}^{(1)}(\alpha) - \frac{3}{8}\alpha b_{3/2}^{(0)}(\alpha), \end{aligned} \quad (2.74) $$

where $\alpha = a_1/a_2$ and the Laplace coefficients are the coefficients of the Laurent expansion (e.g. Laskar, 2005)

$$ (1 - \alpha z)^{-s} (1 - \alpha z^{-1})^{-s} = \frac{1}{2} \sum_{k \in \mathbb{Z}} b_s^{(k)}(\alpha) z^k. \quad (2.75) $$

The variables $X_j$ and $\bar{X}_j$ are not pairs of canonical variables and the equations of motion are not given by Eqs. (2.18). Nevertheless, $(x_j; -i\bar{x}_j)$ is a pair of canonical variables and we can obtain the equations of motion for $X_j$ from those for $x_j$. We have

$$ \dot{X}_j = -\frac{2i}{\Lambda_j} \frac{\partial \mathcal{H}_{\text{LL}}}{\partial \bar{X}_j}, \quad (2.76) $$

and the complete equations of motion take the form

$$ \left(\frac{\dot{X}_1}{\dot{X}_2}\right) = \mathcal{A}\left(\frac{X_1}{X_2}\right) = -\frac{2in_2}{m_0} \begin{pmatrix} m_2\alpha^{-1/2}B_1 & m_2\alpha^{-1/2}B_2 \\ m_1B_2 & m_1B_1 \end{pmatrix} \begin{pmatrix} X_1 \\ X_2 \end{pmatrix}. \quad (2.77) $$

This is a constant linear differential system, whose eigenvalues are the roots of the characteristic polynomial

$$ \det (\mathcal{A} - \lambda \mathbb{I}) = \lambda^2 + \frac{2i B_1 n_2 (m_1 + \alpha^{-1/2} m_2)}{m_0} \lambda + \frac{4 (B_2 - B_1) (B_2 + B_1) m_1 m_2 n_2^2}{m_0^2 \alpha^{1/2}}. \quad (2.78) $$

<footnote n="11">¹¹See the paragraph before Eq. (2.69) in Sect. 2.2.2, for an explanation.</footnote>
For any reasonable choice of masses and semimajor axes, the roots of Eq. (2.78) are pure imaginary. This means that all the trajectories of an averaged, planar, non-resonant two-planet system at second order in eccentricity are such that the semimajor axes are fixed while the eccentricity vectors $e_j e^{i\omega_j}$ behave quasiperiodically (Sect. 2.4.1). Their fundamental frequencies are given by the eigenvalues. For the planetary system Sun-Earth-Jupiter, with the current masses and semimajor axes and assumed planar, a numerical evaluation gives, in arcseconds per year

$$ \lambda_1 = i 7.0758''/\text{yr}, \quad \text{and} \quad \lambda_2 = i 0.009198''/\text{yr}. \qquad (2.79) $$

## 2.4 The differential system $\dot{X} = F(X)$

We review in this section some methods to study differential systems of the form $\dot{X} = F(X)$, where $F: \mathbb{R}^n \to \mathbb{R}^n$ is a differentiable vector field.

### 2.4.1 Fixed points, linearization and eigenvalues

Unless $F$ has a very particular form, the solutions of the Cauchy system

$$ \dot{X} = F(X), \qquad (2.80) $$

where $X(t = 0)$ is given, cannot be written analytically. Nevertheless, if the equation $F(X) = 0$ has solutions on the set of definition of $F$, then these solutions are called fixed points (or equilibria), and in their vicinity, the system (2.80) reduces to a linear system. Let $X_0 \in \mathbb{R}^n$ be such that $F(X_0) = 0$, and we write $X = X_0 + \delta X$. A first-order Taylor expansion in Eq. (2.80) gives

$$ \frac{d\delta X}{dt} = \nabla F(X_0)\delta X + \mathcal{O}(|\delta X|^2), \qquad (2.81) $$

and in the vicinity of $X_0$, the Cauchy system (2.80) is studied by diagonalizing (or reducing to a Jordan form) the matrix $\nabla F(X_0)$. Let $(\lambda_j)_{1 \le j \le n}$ be the eigenvalues of $\nabla F(X_0)$. We call, the linear system, the system (2.81) without $\mathcal{O}(|\delta X|^2)$. If the matrix $\nabla F(X_0)$ is diagonalizable, then the following holds true<fnref n="12" />

* If $\forall j \le n \ \Re \lambda_j \le 0$, then for any initial condition, the flow of the linear system is bounded and the system is stable.

* If $\forall j \le n \ \Re \lambda_j < 0$, then for all initial conditions, the linear system converges towards $X_0$ at exponential speed. The system is asymptotically stable.

* If $\forall j \le n \ \Re \lambda_j = 0$, then all initial conditions lead to a quasiperiodic orbit. If furthermore $\forall i, j \le n \ \Im \lambda_i / \Im \lambda_j \in \mathbb{Q}$, then the trajectories are periodic.

* If $\exists j \le n \ \Re \lambda_j > 0$, then for all initial conditions, but those in a set of measure zero, the linear system diverges away from the fixed point $X_0$ at exponential speed, and its flow is unbounded. The system is unstable.

<footnote n="12">¹²And is easily proven by diagonalizing the system.</footnote>
The condition that $\nabla F(X_0)$ be diagonalizable is necessary. Consider for instance the case where this matrix is nilpotent and non-zero. Then all its eigenvalues are zero but the trajectories are unbounded for almost all initial conditions. We use the results of this section in Sect. 6.2.1 to study the stability, with tides, of a pair of co-orbital planets. We use them again in Sect. 7.2.1 to study the stability, with tides, of the resonance chain $p : p : p + 1$.

## 2.4.2 Eigenvalues of a perturbed matrix

The behaviour of a differential system in the vicinity of its equilibria is described by the knowledge of its eigenvalues. We consider in this subsection the case where the matrix of the linearized differential system takes the form

$$ \mathcal{M} = \mathcal{M}_0 + \varepsilon \mathcal{M}_1, \qquad (2.82) $$

where $\varepsilon$ is a small quantity with respect with unity, written to emphasize the fact that $\mathcal{M}$ is $\varepsilon$-close to $\mathcal{M}_0$. The matrix $\mathcal{M}$ is too complicated for its eigenvalues to be computed directly in an analytical manner, but the perturbative approach we present here, briefly mentioned by Laskar et al. (2012), gives the eigenvalues and eigenvectors very easily. Assume that we know a diagonal basis for $\mathcal{M}_0$

$$ \mathcal{D}_0 = P_0^{-1} \mathcal{M}_0 P_0 = \operatorname{diag}(\lambda_i), \qquad (2.83) $$

where the columns of $P_0$ are the eigenvectors of $\mathcal{M}_0$ and the $\lambda_i$ its eigenvalues, which are not assumed to be of multiplicity one but which are assumed to be sorted by value, that is, equal eigenvalues are consecutive. This does not restrict the generality, as any permutation can be applied on the columns of $P_0$ to achieve that. We now define

$$ \mathcal{Q}_1 = P_0^{-1} \mathcal{M}_1 P_0. \qquad (2.84) $$

Let $P$ denote the matrix of the eigenvectors of $\mathcal{D}_0 + \varepsilon \mathcal{Q}_1$. Since $\mathcal{D}_0 + \varepsilon \mathcal{Q}_1$ is near diagonal, we write

$$ P = \mathbb{I}_n + \varepsilon P_1 + \mathcal{O}(\varepsilon^2). \qquad (2.85) $$

We have

$$ P^{-1} (\mathcal{D}_0 + \varepsilon \mathcal{Q}_1) P = \mathcal{D}_0 + \varepsilon (\mathcal{Q}_1 + [\mathcal{D}_0, P_1]) + \mathcal{O}(\varepsilon^2), \qquad (2.86) $$

where $[\mathcal{D}_0, P_1] = \mathcal{D}_0 P_1 - P_1 \mathcal{D}_0$. Denoting $\mathcal{D}_1 = \operatorname{diag}(q_{i,i})$ the diagonal matrix composed of the diagonal terms of $\mathcal{Q}_1$, we obtain the cohomological equation

$$ \mathcal{Q}_1 + [\mathcal{D}_0, P_1] = \mathcal{D}_1. \qquad (2.87) $$

The solution of the cohomological equation is

$$ p_{i,j} = \begin{cases} \frac{q_{i,j}}{\lambda_j - \lambda_i} & \text{if } |\lambda_i - \lambda_j| \gg \varepsilon, \\ 0 & \text{else,} \end{cases} \qquad (2.88) $$

where $\varepsilon \mathcal{Q}_1 = (q_{i,j})_{1 \le i,j \le n}$ and $\varepsilon P_1 = (p_{i,j})_{1 \le i,j \le n}$. The matrix $\mathcal{M}_0 + \varepsilon \mathcal{M}_1$ is now block diagonal, that is

$$ P^{-1} P_0^{-1} (\mathcal{M}_0 + \varepsilon \mathcal{M}_1) P_0 P = \operatorname{diag}(\mathcal{B}_0^1 + \varepsilon \mathcal{B}_1^1, \dots, \mathcal{B}_0^r + \varepsilon \mathcal{B}_1^r), \quad r \le n \qquad (2.89) $$
where the principal matrix of each block is an homothety

$$ \forall i \le r \quad \exists k \le n \quad \mathcal{B}_0^i = \lambda_k \mathbb{I}_{m(k)}, \qquad (2.90) $$

and $m(k)$ denotes the multiplicity of $\lambda_k$ and thus the size of the block. The computation of the eigenvalues of $M_0 + \varepsilon M_1$ is reduced to the computation of the eigenvalues of the blocks $\mathcal{B}_0^i + \varepsilon \mathcal{B}_1^i$ which are hopefully all of small size and whose eigenvalues are then analytically easily found.

If the eigenvalues of $M_0$ all have multiplicity one and we are only interested in computing the eigenvalues of $M$, then the cohomological Eq. (2.87) does not need to be solved. At first order in $\varepsilon$, the eigenvalues of $M$ can be read on the diagonal of $D_0 + \varepsilon Q_1$. This perturbative method was exposed in appendix E of Couturier et al. (2021). We use it in Sects. 6.2.1 & 6.4.3 to compute the eigenvalues of the differential system describing the motion of two co-orbital planets perturbed by tidal dissipation, in the vicinity of the Lagrangian equilateral equilibria. We also use it in Sect. 3.3.3.

### 2.4.3 Case where $F$ is Hamiltonian

In the case where the vector field $F$ derives from a Hamiltonian function, the differential system linearized in the vicinity of the equilibria has a particular form and this gives constraints on the eigenvalues. Assume that there exists a real Hamiltonian $\mathcal{H}: \mathbb{R}^{2n} \to \mathbb{R}$ such that $F = \mathcal{J}\nabla\mathcal{H}$. In the neighbourhood of a fixed point $X_0$ of the differential system, we note $X = X_0 + \delta X$ and the linearized system reads

$$ \frac{d\delta X}{dt} = \mathcal{J}\nabla^2\mathcal{H}(X_0)\delta X = \mathcal{J}\mathcal{S}\delta X, \qquad (2.91) $$

where $\mathcal{S} = \nabla^2\mathcal{H}(X_0)$ is the Hessian of $\mathcal{H}$ evaluated at $X_0$ (Eq. (2.5)). We now prove that<fnref n="13" /> (Meyer and Hall, 2009, Sect. 3.3)

* If $\lambda$ is eigenvalue of $\mathcal{J}\mathcal{S}$, then so are $-\lambda$, $\bar{\lambda}$ and $-\bar{\lambda}$.

A consequence of this result is that, in the complex plane, the eigenvalues of a linear and real Hamiltonian system are the vertices of rectangles centered at the origin and whose edges are parallel to the axes. Another consequence, in virtue of Sect. 2.4.1, is that a linear and real Hamiltonian system is either unstable or quasiperiodic. In the case $n = 1$, we can also deduce that the two eigenvalues of the system are either $\lambda = \pm\lambda$, or $\lambda = \pm i\lambda$, where $\lambda$ is some positive real number.

We first prove that the matrix $\mathcal{A} = \mathcal{J}\mathcal{S}$ verifies<fnref n="14" /> $\mathcal{A} = \mathcal{J}^t\mathcal{A}\mathcal{J}$. Indeed we have

$$ \mathcal{J}^t\mathcal{A}\mathcal{J} = \mathcal{J}^t(\mathcal{J}\mathcal{S})\mathcal{J} = \mathcal{J}\mathcal{S}^t\mathcal{J}\mathcal{J} = \mathcal{J}\mathcal{S} = \mathcal{A}, \qquad (2.92) $$

where we used the fact that $\mathcal{S}$ is symmetric for the second equality and $^t\mathcal{J}\mathcal{J} = \mathbb{I}$ for the third. Because $\mathcal{H}$ is real, $\mathcal{A}$ is a real matrix and its characteristic polynomial $P(\lambda) = \det(\mathcal{A} - \lambda\mathbb{I})$ is also real. It is thus clear that if $\lambda$ is eigenvalue, so is $\bar{\lambda}$. Proving that $-\lambda$ is eigenvalue

<footnote n="13">¹³The upper bar denotes the complex conjugated.</footnote>

<footnote n="14">¹⁴Such matrices are called Hamiltonian matrices.</footnote>
if $\lambda$ is, will also prove that $-\bar{\lambda}$ is eigenvalue if $\bar{\lambda}$ is. To do that we can show that $\det(\mathcal{A} - \lambda\mathbb{I}) = \det(\mathcal{A} + \lambda\mathbb{I})$. Using $\det \mathbb{J} = 1$, we have

$$ \det (\mathcal{A} - \lambda \mathbb{I}) = \det (\mathbb{J}^t \mathcal{A} \mathbb{J} - \lambda \mathbb{I}) = \det (\mathbb{J}^t \mathcal{A} \mathbb{J} + \lambda \mathbb{J} \mathbb{J}) \qquad (2.93) \\ = \det \mathbb{J} \det (\mathbb{I}^t \mathcal{A} + \lambda \mathbb{I}) \det \mathbb{J} = \det (\mathcal{A} + \lambda \mathbb{I}), $$

which ends the demonstration.

## Protection against destabilization

Another consequence of the previous result is that one degree of freedom Hamiltonians are protected against destabilization by a Hamiltonian perturbation. Assume that a dynamical system is governed by a one degree of freedom Hamiltonian $\mathcal{H}_0(p, q)$, and that $X_0 = (p_0, q_0)$ is a stable equilibrium of that Hamiltonian. According to what was just proven, the eigenvalues in the vicinity of $X_0$ are $\lambda_0 = \pm i\lambda_0$, where we suppose $\lambda_0 > 0$ to prevent trivial dynamics. Assume now that this system is perturbed by a tiny Hamiltonian perturbation

$$ \mathcal{H}(p, q) = \mathcal{H}_0(p, q) + \varepsilon \mathcal{H}_1(p, q), \qquad (2.94) $$

where $\varepsilon \ll 1$ is written to emphasize the smallness of the perturbation. The fixed point is displaced at the position $X = X_0 + \varepsilon X_1$ by the perturbation and the eigenvalues, now given by the spectrum of $\mathbb{J}\nabla^2\mathcal{H}(X)$, are also perturbed to their new values<fnref n="15" />

$$ \lambda^{(1)} = i\lambda_0 + \varepsilon\lambda^{(1)}, \quad \text{and} \quad \lambda^{(2)} = -i\lambda_0 + \varepsilon\lambda^{(2)}, \qquad (2.95) $$

where $\lambda^{(1)}, \lambda^{(2)} \in \mathbb{C}$ are of the same order of magnitude than $\lambda_0$. Since the new differential system is still Hamiltonian, there exists $\lambda \in \mathbb{R}$ such that either $(\lambda^{(1)}, \lambda^{(2)}) = (i\lambda, -i\lambda)$ or $(\lambda^{(1)}, \lambda^{(2)}) = (\lambda, -\lambda)$. Given the small size of $\varepsilon$, the latter possibility is excluded, and we must have $(\lambda^{(1)}, \lambda^{(2)}) = (i\lambda, -i\lambda)$. The new equilibrium is still stable, and one degree of freedom Hamiltonian are protected against destabilization by a Hamiltonian perturbation.

We can use this result to prove that the perturbations due to general relativity (second order in the Post-Newtonian expansion) cannot destabilize a co-orbital system in the same way that tides do. See Sect. 3.3.1 for details.

<footnote n="15">¹⁵Which are conveniently computed in practice using the results of Sect. 2.4.2</footnote>
# The co-orbital motion

Most results of this chapter were first published in Robutel and Pousse (2013) and are recalled in Sect. 2 of Couturier et al. (2021). Section 3.3 is original, while the end of Sect. 3.2.2 was first published in Couturier et al. (2021).

## 3.1 Hamiltonian of two planar co-orbital planets

In this chapter, we consider the point-mass, planar, three-body problem, where a star of mass $m_0$ is orbited by two planets of masses $m_1$ and $m_2$. We introduce the small parameter $\iota = (m_1 + m_2)/m_0$ and 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.

### 3.1.1 Expansion in the neighbourhood of the resonance

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$, we use Poincaré cartesian heliocentric canonical coordinates $(\lambda_j, \tilde{x}_j; \Lambda_j, x_j)$, defined as (Eq. (2.49), Sect. 2.1.4)

$$ \Lambda_j = \beta_j \sqrt{\mu_j a_j}, \quad x_j = \sqrt{\Lambda_j} \sqrt{1 - \sqrt{1 - e_j^2} \exp(i\varpi_j)}, \quad \tilde{x}_j = -i\bar{x}_j. \quad (3.1) $$

With this set of coordinates, the Hamiltonian derives from the symplectic form

$$ \sum_{j \in \{1, 2\}} (d\lambda_j \wedge d\Lambda_j + d\tilde{x}_j \wedge dx_j), \quad (3.2) $$

and according to Eq. (2.54), we write (Robutel and Pousse, 2013)

$$ H = H_K(\Lambda_1, \Lambda_2) + \iota H_P(\Lambda_1, \Lambda_2, \lambda_1, \lambda_2, x_1, x_2, \tilde{x}_1, \tilde{x}_2), \quad (3.3) $$

where the Keplerian part $H_K$, due to star-planet interactions, reads (Eq. (2.55))

$$ H_K(\Lambda_1, \Lambda_2) = - \sum_{j \in \{1, 2\}} \frac{\beta_j^3 \mu_j^2}{2\Lambda_j^2}, \quad (3.4) $$
while the perturbation $\iota H_P$ is due to planet-planet interactions. If the system is at the co-orbital resonance 1 : 1, then the nominal mean motions of both planets comply with $n_1 = n_2 := \eta$. In that case, the semi-major axes $a_j$ of the planets are close to their nominal value $\bar{a}$, defined as

$$ \bar{a} = \mu_0^{1/3} \eta^{-2/3} \quad \text{and} \quad \mu_0 = \mathcal{G} m_0. \qquad (3.5) $$

This also means that the circular angular momenta $\Lambda_j$ are close to their nominal value $\Lambda_j^*$ given by

$$ \Lambda_j^* = m_j \sqrt{\mu_0 \bar{a}}. \qquad (3.6) $$

At the exact mean motion resonance (also called Keplerian resonance), when $n_1 = n_2 = \eta$, the semi-major axes are not equal, but rather verify $a_j = \bar{a}(1 + \mathcal{O}(\iota))$. Considering that the exact co-orbital resonance occurs at $a_1 = a_2 = \bar{a}$ hence generates an error of size $\mathcal{O}(\iota)$, as does substituting $m_j$ for $\beta_j$ or $\mu_0$ for $\mu_j$, in Eqs. (3.5) and (3.6). However, Niederman et al. (2020) have shown that the planets are still in the co-orbital resonance if<fnref n="1" />

$$ \frac{\Lambda_j}{\Lambda_j^*} = 1 + \mathcal{O}(\iota^\xi), \quad \text{where} \quad \frac{1}{3} < \xi \le \frac{1}{2}. \qquad (3.7) $$

In other words, the resonance has a width at least $\mathcal{O}(\iota^{1/2})$ and errors of size $\mathcal{O}(\iota)$ have no consequences. Since we are interested in a study of the dynamics in the vicinity of the resonance, we expand the Hamiltonian in the neighbourhood of $\Lambda_j = \Lambda_j^*$, assuming that $\Lambda_j^*$ is the value of $\Lambda_j$ at the exact resonance.

In the co-orbital resonance, the librating angle $\lambda_1 - \lambda_2$ plays an important role and in order for it to appear explicitly, we perform the canonical transformation (see Eq. (2.31))

$$ \Phi_1(\Lambda_1, \Lambda_2, \lambda_1, \lambda_2) = (Z, Z_2, \phi, \phi_2) = (\Lambda_1 - \Lambda_1^*, \Lambda_1 + \Lambda_2 - \Lambda_1^* - \Lambda_2^*, \lambda_1 - \lambda_2, \lambda_2). \qquad (3.8) $$

At second order in $(Z, Z_2)$, the Keplerian part of the Hamiltonian reads

$$ \hat{H}_K(Z, Z_2) = H_K(\Lambda_1, \Lambda_2) = \eta Z_2 - \frac{3}{2}\eta \left( \frac{Z^2}{\Lambda_1^*} + \frac{(Z_2-Z)^2}{\Lambda_2^*} \right) + \hat{R}_K^{(2)}, \qquad (3.9) $$

where $\hat{H}_K = H_K \circ \Phi_1^{-1}$ and a constant term $-\eta(\Lambda_1^* + \Lambda_2^*)/2$, corresponding to the order 0, has been removed. The remainder $\hat{R}_K^{(2)}$ has a size relative to the expansion $\hat{R}_K^{(2)}/\hat{H}_K = \mathcal{O}(\iota^{2\xi})$, or equivalently<fnref n="2" />

$$ \hat{R}_K^{(2)} = \mathcal{O}(\iota^{3\xi}) H_K. \qquad (3.10) $$

Since the perturbation is small with respect to the Keplerian part, it is not necessary to expand it up to the second order in the vicinity of $\Lambda_j^*$. If we limit ourselves to the zeroth order, that is

$$ H_P(\Lambda_1, \Lambda_2, \lambda_1, \lambda_2, x_1, x_2, \tilde{x}_1, \tilde{x}_2) = H_P(\Lambda_1^*, \Lambda_2^*, \lambda_1, \lambda_2, x_1, x_2, \tilde{x}_1, \tilde{x}_2) + R_P^{(0)}, \qquad (3.11) $$

then the size of the remainder relative to the expansion is $R_P^{(0)}/H_P = \mathcal{O}(\iota^{\xi})$, that is

$$ R_P^{(0)} = \mathcal{O}(\iota^{1+\xi}) H_K. \qquad (3.12) $$

<footnote n="1">¹This domain for $\xi$ corresponds to the horseshoe-shaped orbits (Sect. 3.2.1).</footnote>

<footnote n="2">²The removal of the constant $-\eta(\Lambda_1^* + \Lambda_2^*)/2$ from $\hat{H}_K$ yields $\hat{H}_K/H_K = \mathcal{O}(\iota^{\xi})$.</footnote>
Combining the estimates (3.10) and (3.12), we obtain

$$ \frac{R_P^{(0)}}{\hat{R}_K^{(2)}} = \mathcal{O}(\iota^{1-2\xi}). \qquad (3.13) $$

We can verify that truncating the expansion (3.11) at order zero in the vicinity of the Keplerian resonance is enough. When $\xi \le 1/2$, Eq. (3.13) shows that $R_P^{(0)}$ is negligible with respect to $\hat{R}_K^{(2)}$, whereas when $\xi > 1/2$ (this corresponds to tadpole orbits, see Sect. 3.2.1), the estimates (3.10) and (3.12) show that both remainders quickly tend to 0 with increasing $\xi$. In both cases, the order 0 is enough and we simply evaluate the perturbation $H_P$ at $\Lambda_j = \Lambda_j^*$.

Lastly, for the expansion (3.9) to be valid, it is important that $\hat{R}_K^{(2)}$ is negligible with respect to $H_P$. Combining $R_P^{(0)}/H_P = \mathcal{O}(\iota^{\xi})$ with Eq. (3.13), we can write $\hat{R}_K^{(2)}/H_P = \mathcal{O}(\iota^{3\xi-1})$, and we must have $\xi > 1/3$. This is the lower bound of Niederman et al. (2020).

### 3.1.2 Averaged Hamiltonian

The timescales of evolution can be established by uncoupling the action variables in the Keplerian Hamiltonian. We perform the canonical transformation (Eq. (2.31))

$$ \Phi_2(Z, Z_2, \phi, \phi_2) = (I, I_2, \xi, \xi_2) = (Z - \delta Z_2, Z_2, \phi, \delta\phi + \phi_2), \quad \text{with} \quad \delta = \frac{m_1}{m_1 + m_2}, \qquad (3.14) $$

where $x_j$ and $\tilde{x}_j$ are unchanged. The Keplerian part of the Hamiltonian takes the form

$$ \check{H}_K(I, I_2) = \hat{H}_K(Z, Z_2) = \eta I_2 - \frac{3}{2}\eta \frac{I_2^2}{\Lambda_1^* + \Lambda_2^*} - \frac{3}{2}\eta \frac{\Lambda_1^* + \Lambda_2^*}{\Lambda_1^* \Lambda_2^*} I^2, \qquad (3.15) $$

where $\check{H}_K = \hat{H}_K \circ \Phi_2^{-1}$. This decoupling emphasizes the timescales involved in the dynamics. The Hamilton Eqs. (2.18) give

$$ \dot{\xi}_2 = \frac{\partial \check{H}_K}{\partial I_2} = \eta (1 + \mathcal{O}(\iota^\xi)) \quad \text{and} \quad \dot{\xi} = \frac{\partial \check{H}_K}{\partial I} = \eta \mathcal{O}(\iota^\xi), \qquad (3.16) $$

and so, the angle $\xi_2$ evolves much more quickly than $\xi = \lambda_1 - \lambda_2$. As we will show in Sect. 3.2.2, the eccentricities $x_j$ evolves on timescales much longer than the timescale of evolution of $\xi$. This means that the co-orbital dynamics takes place on three distinct timescales. The fast one, of frequency $\eta \mathcal{O}(1)$ is associated with $\xi_2$ (or either of the $\lambda_j$), the semi-fast one, of frequency $\eta \mathcal{O}(\iota^\xi)$, is associated with the libration of $\xi$, while the secular (or slow) one has a frequency $\eta \mathcal{O}(\iota)$ and is associated with the $x_j$. This will be further detailed in Sects. 3.2.1 & 3.2.2.

Since the angle $\xi_2$ evolves at the orbital frequency, it is natural to average the Hamiltonian over it; that is, we perform a first-order Lie serie expansion (Sect. 2.2.2) and write$^3$

$$ \check{H}(I, I_2, \xi, x_1, x_2, \tilde{x}_1, \tilde{x}_2) = \check{H}_K(I, I_2) + \frac{1}{2\pi} \int_0^{2\pi} \iota \check{H}_P(\xi, \xi_2, x_1, x_2, \tilde{x}_1, \tilde{x}_2) d\xi_2 + \check{R}^{(2)}, \qquad (3.17) $$

<footnote>$^3$The Lie expansion performs on the Hamiltonian a canonical transformation $\iota$-close to the identity, but for convenience, we do not change the name of the variables.</footnote>
where $\tilde{H}_P = H_P \circ \Phi_1^{-1} \circ \Phi_2^{-1}$ and $\tilde{R}^{(2)}$ is the remainder of the Lie expansion. Requesting that $\tilde{R}^{(2)} \ll \tilde{H}_P$ yields once again $\xi > 1/3$ (Niederman et al., 2020). After the averaging process, the Hamiltonian no longer depends on $\xi_2$ and its conjugated momentum $I_2$ is a first integrale. That is, the total circular angular momentum $\Lambda_1 + \Lambda_2$ is conserved in the averaged problem. In the co-orbital resonance 1 : 1, this quantity corresponds, for other resonances, to what various authors call the scaling parameter (Michtchenko et al., 2008; Delisle, 2017; Petit et al., 2020). In the averaged problem, the conservation of $\Lambda_1 + \Lambda_2$ allows us to lose one degree of freedom, and we are left with the three degrees of freedom $(\xi, \tilde{x}_1, \tilde{x}_2; I, x_1, x_2)$.

While the total circular angular momentum $\Lambda_1 + \Lambda_2$ is conserved in the averaged problem, it is well known that the total angular momentum $\sum_j \Lambda_j - |x_j|^2$ is also conserved$^4$, even in the complete (non-averaged) problem. Hence, there exists a transformation allowing one more degree of freedom to be lost, and the averaged planar co-orbital problem has two degrees of freedom. More generally, an averaged, resonant, two-planet planar system has two degrees of freedom$^5$. The transformation to two degrees of freedom is performed by Giuppone et al. (2010) in their numerical study of the co-orbital motion. However, the subsequent set of coordinates is singular at zero eccentricity, which is not adapted for an analytical work. Hence we stick here with three degrees of freedom.

In order to work with dimensionless quantities, we normalise the action variables by

$$J = \frac{I}{m\bar{a}^2\eta}, \quad J_2 = \frac{I_2}{m\bar{a}^2\eta}, \quad X_j = \sqrt{\frac{2}{m_j\bar{a}^2\eta}} x_j, \quad \bar{X}_j = i\sqrt{\frac{2}{m_j\bar{a}^2\eta}} \tilde{x}_j, \quad (3.18)$$

where $m = \sqrt{m_1 m_2}$ and the upper bar denotes the complex conjugated. These new eccentricity variables comply with

$$X_j = e_j \exp(i\varpi_j) \left(1 + \mathcal{O}(\iota^5) + \mathcal{O}(e_j^2)\right). \quad (3.19)$$

The order of magnitude of $J$ and $J_2$ is $\mathcal{O}(\iota^5)$ and these variables are linked to the semi-major axes $a_j$ by

$$a_j = \bar{a} \mathcal{R}_j^2, \quad (3.20)$$

with

$$\mathcal{R}_j = 1 + f_j, \quad f_1 = \frac{m}{m_1 + m_2} J_2 + \frac{m}{m_1} J, \quad f_2 = \frac{m}{m_1 + m_2} J_2 - \frac{m}{m_2} J. \quad (3.21)$$

The transformations $\Phi_1$ and $\Phi_2$ are canonical (Sect. 2.1.3), but the normalization (3.18) is not. We can nevertheless have equations of motion close to their canonical form (2.18) by rescaling the Hamiltonian by

$$\mathcal{H} = \frac{\tilde{H}}{m\bar{a}^2\eta}, \quad (3.22)$$

and we obtain, for the equations of motion (see Eq. (2.76))

$$\dot{J} = -\frac{\partial \mathcal{H}}{\partial \xi}, \quad \dot{J}_2 = -\frac{\partial \mathcal{H}}{\partial \xi_2}, \quad \dot{\xi} = \frac{\partial \mathcal{H}}{\partial J}, \quad \dot{\xi}_2 = \frac{\partial \mathcal{H}}{\partial J_2}, \quad \dot{X}_j = -2i \frac{m}{m_j} \frac{\partial \mathcal{H}}{\partial \bar{X}_j}, \quad \dot{\bar{X}}_j = \bar{X}_j. \quad (3.23)$$

<footnote>$^4$The total AMD $\sum_j |x_j|^2$ (Laskar, 2000) is thus conserved in the averaged 1 : 1 resonance.</footnote>

<footnote>$^5$Henrard and Lemaitre (1983) even showed that at first order in eccentricity, the system has only one degree of freedom.</footnote>
With this last transformation, the Keplerian part of the Hamiltonian reads

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

where the constant terms depending only on the first integral $J_2$ have been ignored. In Chap. 6, due to tidal dissipation, $J_2$ is no longer a first integral and $\mathcal{H}_K$ includes more terms (see Eq. (6.12)).

### 3.1.3 Expansion of the perturbative part

In the set of canonical coordinates $(\tilde{\mathbf{r}}_j; \mathbf{r}_j)$, where $\mathbf{r}_j$ and $\tilde{\mathbf{r}}_j$ are respectively the heliocentric positions and barycentric momenta, the perturbation to the Hamiltonian takes the form given by Eq. (2.53)

$$ \iota H_P(\mathbf{r}_j, \tilde{\mathbf{r}}_j) = \frac{\tilde{\mathbf{r}}_1 \cdot \tilde{\mathbf{r}}_2}{m_0} - \frac{Gm_1m_2}{|\mathbf{r}_1 - \mathbf{r}_2|}. \qquad (3.25) $$

We write<fnref n="6" />

$$ |\mathbf{r}_1 - \mathbf{r}_2|^2 = r_1^2 + r_2^2 - 2r_1r_2 \cos(\mathbf{r}_1, \mathbf{r}_2) = \bar{a}^2 \Delta^2 (1 + \mathcal{O}(e_j^2)), \qquad (3.26) $$

where $\Delta = \sqrt{2 - 2\cos\xi}$, and we expand $H_P$ in power series of the eccentricities. The expansion takes the form (Robutel and Pousse, 2013; Couturier et al., 2021)

$$ \iota H_P = \sum_{n=0}^{+\infty} \mathcal{H}^{(2n)}, \quad \text{with} \quad \mathcal{H}^{(2n)} = \sum_{|\mathbf{p}|=2n} \Psi_{\mathbf{p}}(\Delta^{-1}, e^{i\xi}, e^{-i\xi}) X_1^{p_1} X_2^{p_2} \bar{X}_1^{p_3} \bar{X}_2^{p_4}, \qquad (3.27) $$

where the norm of the tuple $\mathbf{p} = (p_1, p_2, p_3, p_4) \in \mathbb{N}^4$ is defined by $|\mathbf{p}| = p_1 + p_2 + p_3 + p_4$ and $\Psi_{\mathbf{p}}$ is polynomial in its arguments. The d'Alembert rule (Eq. (2.59))

$$ p_1 + p_2 = p_3 + p_4, \qquad (3.28) $$

equivalent to the conservation of the total angular momentum, justifies that the expansion has no terms of odd order in eccentricity. In practice, the expansion (3.27) can be explicitly computed from the results of Laskar and Robutel (1995). We apply the transformations $\Phi_1$ and $\Phi_2$ to $H_P$, as well as the normalisation (3.18) and the rescaling (3.22). We use the algebraic manipulator *Trip* (Gastineau and Laskar, 2011), and up to the fourth order in eccentricity, we obtain (Couturier et al., 2021)

$$ \begin{aligned} \mathcal{H}^{(0)} &= \frac{m}{m_0} \eta (\cos \xi - \Delta^{-1}), \\ \mathcal{H}^{(2)} &= \frac{1}{2} \frac{m}{m_0} \eta \{A_h (X_1 \bar{X}_1 + X_2 \bar{X}_2) + B_h X_1 \bar{X}_2 + \bar{B}_h \bar{X}_1 X_2\}, \\ \mathcal{H}^{(4)} &= \frac{1}{4} \frac{m}{m_0} \eta \{D_h (X_1^2 \bar{X}_1^2 + X_2^2 \bar{X}_2^2) + E_h X_1^2 \bar{X}_2^2 + \bar{E}_h X_2^2 \bar{X}_1^2 \\ &\quad + F_h (X_1 X_2 \bar{X}_1^2 + \bar{X}_1 X_2 \bar{X}_2^2) + \bar{F}_h (\bar{X}_1 \bar{X}_2 X_1^2 + X_1 X_2 \bar{X}_2^2) + G_h X_1 X_2 \bar{X}_1 \bar{X}_2\}, \end{aligned} \qquad (3.29) $$

<footnote n="6">⁶This is because we expand at the zeroth order in $\Lambda_j$, see Sect. 3.1.1.</footnote>
where the coefficients of the second order are

$$ 
\begin{aligned}
A_h &= \frac{5 \cos 2\xi - 13 + 8 \cos \xi}{4\Delta^5} - \cos \xi \quad \text{and} \\
B_h &= e^{-2i\xi} - \frac{e^{-3i\xi} + 16e^{-2i\xi} - 26e^{-i\xi} + 9e^{i\xi}}{8\Delta^5}.
\end{aligned} 
\quad (3.30) 
$$

The coefficients $D_h$, $E_h$, $F_h$ and $G_h$ of the fourth order in eccentricity are given in Eqs. (B.1) of appendix B.

## 3.2 Lagrangian equilibria and linearization

### 3.2.1 Circular dynamics

The expansion of the perturbative part of the Hamiltonian has no terms of odd degree in the variables $X_j$ (Eq. (3.29)). As a consequence, the form of $dX_j/dt$ in Eqs. (3.23) implies that the manifold $X_1 = X_2 = 0$ is stable by the flow of the averaged Hamiltonian. In other words, initially circular orbits remain circular. The dynamics on this stable manifold is driven by the one-degree-of-freedom Hamiltonian $\mathcal{H}_K + \mathcal{H}^{(0)}$, whose phase space we plot in Fig. 3.1. This Hamiltonian has three fixed points at $(J = 0, \xi = \pi/3)$, $(J = 0, \xi = 5\pi/3)$ and $(J = 0, \xi = \pi)$, corresponding to the famous Lagrangian equilibria. The two elliptic points $L_4$ and $L_5$, located on energy maximizers, correspond to configurations where the three bodies are at the vertices of an equilateral triangle, while the hyperbolic configuration $L_3$ is a saddle point, where the bodies are aligned. The separatrix emanating from $L_3$

<figure><img src="image_5.png" /><figcaption>Fig. 3.1 — Phase portrait of the Hamiltonian $\mathcal{H}_K + \mathcal{H}^{(0)}$ with $m_1 = 0.001 m_0$ and $m_2 = 0.0004 m_0$. The two elliptic fixed points $L_4$ and $L_5$ are energy maximizers, while the hyperbolic fixed point $L_3$ is a saddle point. Around $L_{4,5}$ are the tadpoles orbits and outside of the separatrix emanating from $L_3$ are the horseshoe-shaped orbits. The Hamiltonian diverges to minus infinity as $\lambda_1 - \lambda_2 = \xi$ tends to zero.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 3.2 — Relative position of $m_1$ with respect to $m_2$ along the trajectories of the Hamiltonian $\mathcal{H}_K + \mathcal{H}^{(0)}$. The trajectories are numerically integrated with a Dormand-Prince method of order 8 using initial conditions $J = 0$ and either $\xi = \pm 45^\circ$ (orange tadpoles) or $\xi = 10^\circ$ (black horseshoe). The curves show the motion of $m_1$ relative to the black dot $m_2$; that is, the $x$ and $y$ axes of this figure are $(x, y) = (a_1 \cos \xi - a_2, a_1 \sin \xi)$, the $a_j$ coming from Eq. (3.20) with $J_2 = 0$. The co-orbital masses are $m_2 = 5 m_1 = 0.001 m_0$. The period of the motion is here $12.11 T$ for the black curve and $11.32 T$ for the orange ones, where $T = 2\pi/\eta$.</figcaption></figure>

splits the phase space in three regions. Around the elliptic fixed points $L_4$ and $L_5$ are orbits commonly known as tadpole, while the orbits encircling $L_3$ are called horseshoe-shaped orbits, due to the shape of the trajectories in the rotating frame (see Fig. 3.2). In our model, the angle $\xi = \lambda_1 - \lambda_2$ is always librating and the exterior of the resonance does not exist. If no expansion in the vicinity of the Keplerian resonance $\Lambda_j = \Lambda_j^*$ is performed, the phase space contains two more hyperbolic fixed points labeled $L_1$ and $L_2$, and the separatrices emanating from them are the boundary of the co-orbital resonance (see Fig. 1 of Robutel and Pousse, 2013). Instead, our phase space features a singular line $\xi = 0$ where the Hamiltonian diverges to $-\infty$, and the model is not valid for very large libration amplitudes.

We show with Eq. (3.16) that the libration frequency of $\xi$ is $\eta \mathcal{O}(\iota^{\xi})$. In the vicinity of $L_{4,5}$, we can give an exact expression by linearizing the differential system. We denote $\Delta\xi = \xi - \pi/3$ and in the neighbourhood of $L_4$ (the dynamics are symmetric in the neighbourhood of $L_5$), the linearized system reads

$$ \frac{d}{dt} \begin{pmatrix} \Delta \xi \\ J \end{pmatrix} = \eta \begin{pmatrix} 0 & -3\iota \frac{m_0}{m} \\ \frac{9m}{4m_0} & 0 \end{pmatrix} \begin{pmatrix} \Delta \xi \\ J \end{pmatrix}, \quad (3.31) $$

and its eigenvalues are $\pm i\nu$ where

$$ \nu = \eta \sqrt{\frac{27\iota}{4}} \quad (3.32) $$

is the libration frequency of $\xi$ in the vicinity of the elliptic Lagrangian equilibria.

For a given trajectory in the phase space, we denote $\xi_{cr}$, respectively $J_{cr}$, the extremal value reached by $\xi$, respectively by $J$, along that trajectory. When the action is $J = J_{cr}$, the angle is $\xi = \pi/3$, whereas $\xi = \xi_{cr}$ implies $J = 0$ (Fig. 3.1). The conservation of
$\mathcal{H}_K + \mathcal{H}^{(0)}$ (Eq. (2.32)) hence gives the relation

$$-\frac{3(m_1 + m_2)m_0}{2m_1m_2}J_{\mathbf{e}_x}^2 - \frac{1}{2} = \cos \xi_{\mathbf{e}_x} - (2 - 2 \cos \xi_{\mathbf{e}_x})^{-1/2}, \quad (3.33)$$

that holds whether the extremum is a minimum or a maximum. The boundary between tadpole and horseshoe-shaped orbits is the separatrix emanating from $L_3$, and for the largest tadpole orbit (or the the smallest horseshoe-shaped orbit), $\xi_{\mathbf{e}_x} = \pi$ yields $J_{\mathbf{e}_x} = \mathcal{O}(\iota^{1/2})$. The boundary between the tadpole and horseshoe domain thus corresponds<fnref n="7" /> to $\xi = 1/2$, which is the upper bound in Eq. (3.7). For $\xi < 1/2$, the system is in the horseshoe domain while it is in the tadpole domain for $\xi \ge 1/2$.

For a horseshoe-shaped orbit with a very large libration amplitude, we call $\xi_{\min}$ the minimal value reached by $\xi$, that is, the minimal angular separation between both co-orbitals. In the following, we compute the value of $\xi_{\min}$ such that the co-orbital planets stay outside of each other Hill's spheres. The closest approach $\mathbf{d}$ is given by

$$\mathbf{d} = \bar{a}\xi_{\min}. \quad (3.34)$$

At the closest approach, we perform a Laurent expansion around $\xi_{\min} = 0$ in the right-hand side of Eq. (3.33), retaining only terms of order $-1$ and $0$ in $\xi_{\min}$. This gives

$$\frac{3}{2}\mathbf{d}\left(1 + \frac{\iota m_0^2}{m_1 m_2} J_{\mathbf{e}_x}^2\right) = 1. \quad (3.35)$$

For a large horseshoe-shaped orbit, $\xi < 1/2$ implies that $1 \ll \iota m_0^2 J_{\mathbf{e}_x}^2 / m_1 m_2$ and the estimate $J_{\mathbf{e}_x} = \mathcal{O}(\iota^\xi)$ yields

$$\mathbf{d} = \mathcal{O}(\iota^{1-2\xi}) \bar{a}. \quad (3.36)$$

Our model is based on a perturbative approach where the planet–planet interactions are assumed to be much smaller than the star–planet interactions. For this approach to make sense, it is necessary that the co-orbital bodies stay well outside of each other Hill's spheres, that is

$$\mathbf{d} \gg \bar{a} \left(\frac{\iota}{3}\right)^{1/3}, \quad (3.37)$$

or equivalently, using Eq. (3.36)

$$\mathcal{O}(\iota^{2/3-2\xi}) \gg 1. \quad (3.38)$$

We retrieve once again the lower bound $1/3 < \xi$ given by Niederman et al. (2020), showing that this lower bound can be obtained in three different manners. We can either request that

* the remainder $\hat{R}_K^{(2)}$ of the expansion (3.9) around the Keplerian resonance is negligible with respect to the perturbative part $H_P$ of the Hamiltonian, or that

* the remainder $\tilde{H}^{(2)}$ of the Lie serie expansion (3.17) is negligible with respect to $\tilde{H}_P$, or that

* the co-orbital planets do not enter each other Hill's spheres.

These three conditions are equivalent to $1/3 < \xi$ in Eq. (3.7).

<footnote n="7">⁷This is due to $J_{\mathbf{e}_x} = \mathcal{O}(\iota^\xi)$.</footnote>
### 3.2.2 Eccentric dynamics

The secular dynamics of the problem is driven by the Hamiltonian $\mathcal{H}_K + \mathcal{H}^{(0)} + \sum_{n \ge 1} \mathcal{H}^{(2n)}$. If we limit the summation to $n = 1$; that is, if we only consider $\mathcal{H} = \mathcal{H}_K + \mathcal{H}^{(0)} + \mathcal{H}^{(2)}$, the equations of motion are linear in $X_j$ (see Eqs. (3.23) & (3.29)), which is convenient to obtain analytical results at low eccentricity. We show in Chap. 6 that tidal dissipation damps the eccentricities to 0, and the analytical insight at low eccentricity allows us to understand most of the dynamics. Hence we only consider $\mathcal{H} = \mathcal{H}_K + \mathcal{H}^{(0)} + \mathcal{H}^{(2)}$ in the rest of this chapter and we will include $\mathcal{H}^{(4)}$ back again in Sect. 6.1.1 when we include tidal dissipation in the analytical model. At this order in eccentricity, the equations of motion read

$$ \begin{pmatrix} \dot{X}_1 \\ \dot{X}_2 \end{pmatrix} = -\frac{i\eta}{m_0} \begin{pmatrix} m_2 A_h(\xi(t)) & m_2 \bar{B}_h(\xi(t)) \\ m_1 B_h(\xi(t)) & m_1 A_h(\xi(t)) \end{pmatrix} \begin{pmatrix} X_1 \\ X_2 \end{pmatrix}, \qquad (3.39) $$

where $(J(t), \xi(t))$ is a solution of $\mathcal{H}_K + \mathcal{H}^{(0)}$, such as those plotted in Fig. 3.2. The system of Eqs. (3.39) is linear and periodic of time and in general, its stability can be studied, at least numerically, with the Floquet theory. Nevertheless, at one of the Lagrangian equilibria, the solutions of $\mathcal{H}_K + \mathcal{H}^{(0)}$ are constant and system (3.39) becomes both linear and time-independent, and is studied by diagonalisation. The hyperbolic point $L_3$ is unstable, and we are not interested in its dynamics. Both elliptic equilibria $L_4$ and $L_5$ are symmetric by the transformation $\xi \mapsto -\xi$ and it is enough to study one of them only. We study the eccentric dynamics at $L_4$ and we set $(J, \xi) = (0, \pi/3)$ in Eq. (3.39); that is, we set (see Eq. (3.30))

$$ A_h = -\frac{27}{8} \quad \text{and} \quad B_h = \frac{27}{8}e^{-i\pi/3}. \qquad (3.40) $$

The matrix

$$ M_0 = -\frac{i\eta}{m_0} \begin{pmatrix} m_2 A_h & m_2 \bar{B}_h \\ m_1 B_h & m_1 A_h \end{pmatrix} \qquad (3.41) $$

of the differential system has eigenvalues $ig_1$ and $ig_2$ where

$$ ig_1 = i\eta \frac{27i}{8}, \quad ig_2 = 0, \qquad (3.42) $$

with associated eigenvectors

$$ V_1 = \begin{pmatrix} m_2 e^{4i\pi/3} \\ m_1 \end{pmatrix} \quad \text{and} \quad V_2 = \begin{pmatrix} e^{i\pi/3} \\ 1 \end{pmatrix}. \qquad (3.43) $$

The neutral eigenmode, collinear to $V_2$, is the so-called Lagrange configuration where

$$ \frac{X_1}{X_2} = \exp\left(i\frac{\pi}{3}\right). \qquad (3.44) $$

Since the corresponding eigenvalue $ig_2$ is zero, this configuration corresponds to a family of fixed points parameterized by $J_2 \propto \Lambda_1 + \Lambda_2$, where the eccentricities verify $e_1 = e_2$ and $\varpi_1 - \varpi_2 = \pi/3$. The precessing eigenmode, collinear to $V_1$, is called anti-Lagrange configuration by Giuppone et al. (2010) and verify

$$ \frac{X_1}{X_2} = \frac{m_2}{m_1} \exp\left(i\frac{4\pi}{3}\right). \qquad (3.45) $$
Along this family of periodic orbits, the eccentricities are constant such that $m_1e_1 = m_2e_2$ and the pericentres precess at the same frequency $\eta g_1$ while maintaining the relation $\varpi_1 - \varpi_2 = 4\pi/3$.

For an arbitrary solution $(J(t), \xi(t))$, the differential system (3.39) has a priori to be studied with the Floquet theory, which generally does not allow for analytical expressions<fnref n="8" />. Nevertheless, since $g_1 \ll \nu$, the system (3.39) can be studied by time-averaging over one period of $\xi$. This is justified mathematically by the existence of a canonical transformation $\Psi$ (Theorem 4.5 of Niederman et al., 2020, second averaging theorem)

$$ \Psi(J, \xi) = (J, \xi), \qquad (3.46) $$

such that the dependency on $\xi$ is pushed to the second order in eccentricity, that is

$$ \mathcal{H}(J, \xi, X_j, \bar{X}_j) = \mathcal{H}(J, \xi, X_j, \bar{X}_j) = \mathcal{H}_K(J) + \mathcal{H}^{(0)}(J) + \mathcal{H}^{(2)}(J, \xi, X_j, \bar{X}_j), \quad (3.47) $$

where $\mathcal{H} = \mathcal{H} \circ \Psi^{-1}$. A first-order Lie serie expansion (Sect. 2.2.2) performed on Eq. (3.47) would then allow the dependency on $\xi$ to be lost, where $\mathcal{H}_K + \mathcal{H}^{(0)}$ plays the role of $\mathcal{H}_0$ in Eq. (2.60).

The coefficients of the time average of $\mathcal{M}_0(t)$ cannot be given analytically, but in the horseshoe domain, the symmetry of the orbits allows some analytical insight. The expression of $B_h$ in Eq. (3.30) shows that $\mathfrak{Im}(B_h(2\pi - \xi)) = -\mathfrak{Im}(B_h(\xi))$, and since the trajectories $(J(t), \xi(t))$ are symmetric around $(0, \pi)$ in the horseshoe domain (see Fig. 3.1), the time average of $\mathfrak{Im} B_h$ vanishes for horseshoe-shaped orbits. As a consequence, the coefficients of the averaged matrix

$$ \underline{\mathcal{M}}_0 = \frac{1}{2\pi} \int_0^{2\pi} \mathcal{M}_0(\xi) d\xi = -\frac{i\eta}{m_0} \begin{pmatrix} m_2 \underline{A}_h & m_2 \underline{B}_h \\ m_1 \underline{B}_h & m_1 \underline{A}_h \end{pmatrix}, \quad (3.48) $$

are pure imaginary. Diagonalising $\underline{\mathcal{M}}_0$ shows that its eigenvectors are such that

$$ \arg \left( \frac{X_1}{X_2} \right) = \varpi_1 - \varpi_2 = 0, \qquad (3.49) $$

for the Lagrange-like configuration, while

$$ \arg \left( \frac{X_1}{X_2} \right) = \varpi_1 - \varpi_2 = \pi, \qquad (3.50) $$

for the anti-Lagrange-like configuration. This corresponds to aligned and anti-aligned pericentres. We show in Chap. 6 that for a system initially close to $L_4$, tidal dissipation induces an exponential increase of the libration amplitude, and the system ends up reaching the horseshoe-shaped orbits. In that case, for a system in the Lagrange configuration (resp. in the anti-Lagrange configuration), the difference of the pericentres $\varpi_1 - \varpi_2$ starts at $\pi/3$ (resp. $4\pi/3$) and decreases until it reaches 0 (resp. $\pi$) when the system reaches the horseshoe-shaped orbits. This is very visible in Fig. 6.4

<footnote n="8">⁸The monodromy matrix cannot be computed analytically.</footnote>
## 3.3 Contribution of general relativity

Up to this point, we have studied the point-mass co-orbital resonance 1 : 1 according to the Newtonian laws. In this section, we maintain the point-mass approximation but we now consider that the three bodies interact according to the laws of general relativity. As long as the nominal semimajor axis $\bar{a}$ is large with respect to the Schwarzschild radius

$$r_s = \frac{2Gm_0}{c^2} \qquad (3.51)$$

of the star, where c is the speed of light, the general relativity only slightly perturbs the motion as described by the laws of Newton, and we will only consider in this section the first Post-Newtonian corrections of the general relativity.

### 3.3.1 Some intuition

In the circular case, the dynamics of the co-orbital resonance in the Newtonian, point-mass approximation are described by the one degree of freedom Hamiltonian (Eqs. (3.24) & (3.29))

$$\mathcal{H}_0(J, \xi) = \mathcal{H}_K(J) + \mathcal{H}^{(0)}(J, \xi). \qquad (3.52)$$

Since the perturbations due to general relativity can be written in Hamiltonian form (Eq. (3.54)), there exists a Hamiltonian $\mathcal{H}_1(J, \xi)$ such that, in the circular case, the Hamiltonian of the co-orbital resonance in the Post-Newtonian expansion reads

$$\mathcal{H}_0(J, \xi) + \varepsilon\mathcal{H}_1(J, \xi), \qquad (3.53)$$

where $\varepsilon = v^2/c^2 \approx r_s/(2\bar{a}) \ll 1$ is written to emphasize the smallness of $\mathcal{H}_1$ with respect to $\mathcal{H}_0$ and $v = \bar{a}\eta$ is the orbital speed of the co-orbitals. The two Lagrangian equilibria of $\mathcal{H}_0$ at $(J, \xi) = (0, \pm\pi/3)$, are elliptic (stable), since the eigenvalues in their vicinity, $\pm i\nu$ (see Eq. (3.32)), are pure imaginary. The eigenvalues of $\mathcal{H}_0 + \varepsilon\mathcal{H}_1$ will be a slight perturbation of $\pm i\nu$, and the Lagrangian equilibria may become unstable under the influence of general relativity if one of the eigenvalues is perturbed to a strictly positive real part (Sect. 2.4.1).

Nevertheless, we proved at the end of Sect. 2.4.3 that one degree of freedom Hamiltonians are protected against destabilization by a Hamiltonian perturbation. This means that the eigenvalues in the vicinity of the Lagrangian equilibria will indeed be perturbed by general relativity, but the perturbations will themselves be pure imaginary, and we know without calculations that the Lagrangian equilibria are stable under the effect of general relativity.

This is at least true at the first post-Newtonian approximation, because at the fifth order in $v/c$ and beyond, the emission of gravitational waves induces a dissipation of the orbital energy and prevents the perturbations due to general relativity to be written in Hamiltonian form. Since the Lagrangian equilibria $L_4$ and $L_5$ are energy maximizers (Fig. 3.1), a dissipation of the total orbital energy should in that case make the Lagrangian equilibria unstable.
### 3.3.2 The Post-Newtonian Hamiltonian

At first order in $\epsilon = v^2/c^2$, the Hamiltonian of two bodies of masses $m_0$ and $m_j$ reads (e.g. Landau and Lifshitz, 1980, page 366)

$$ \mathcal{H}_{RG}(\hat{\mathbf{r}}_j, \mathbf{r}_j) = \frac{\hat{\mathbf{r}}_j^2}{2\beta_j} - \frac{\mathcal{G}m_0m_j}{r_j} - \frac{\hat{\mathbf{r}}_j^4}{8c^2}\left(\frac{1}{m_0^3} + \frac{1}{m_j^3}\right) + \frac{\mathcal{G}^2m_0m_j(m_0+m_j)}{2c^2r_j^2} \\ - \frac{\mathcal{G}}{2c^2r_j}\left\{3\hat{\mathbf{r}}_j^2\left(\frac{m_j}{m_0} + \frac{m_0}{m_j} + \frac{7}{3}\right) + \left(\hat{\mathbf{r}}_j \cdot \frac{\mathbf{r}_j}{r_j}\right)^2\right\}, \qquad (3.54) $$

where $\hat{\mathbf{r}}_j$ and $\mathbf{r}_j$ are respectively the barycentric momentum and heliocentric position, defined in Sect. 2.1.4. The first two terms give the Hamiltonian of the two-body problem in the Newtonian approximation (see Eq. (2.36)), and we remove them in order to only keep the relativistic perturbations. Due to the low mass of the planets with respect to that of the star, we do not take into account the planet–planet non-Newtonian interactions. Neglecting $m_j/m_0$ in front of 1, we write the relativistic perturbations to the Hamiltonian of the co-orbital problem as $\epsilon\mathcal{H}_1 = \mathcal{H}_{RG}^{(1)} + \mathcal{H}_{RG}^{(2)}$, where

$$ \mathcal{H}_{RG}^{(j)} = -\frac{\hat{\mathbf{r}}_j^4}{8c^2m_j^3} + \frac{\mathcal{G}^2m_0^2m_j}{2c^2r_j^2} - \frac{3\mathcal{G}m_0\hat{\mathbf{r}}_j^2}{2m_jc^2r_j}. \qquad (3.55) $$

As in Sect. 3.1.3, we rewrite the Hamiltonian $\mathcal{H}_{RG}^{(j)}$ in the variables $(J, J_2, \xi, \xi_2, X_j, \bar{X}_j)$. We are only interested in the secular (i.e. long-term) evolution of the orbits, and we average over the angle $\xi_2$. Even though the perturbation $H_P$, due to planet–planet Newtonian interactions, was expanded at the zeroth order (Eq. (3.11)) in the variables $\Lambda_j$ around the Keplerian resonance (Eq. (3.6)), we keep exact expressions<fnref n="9" /> in the $\Lambda_j$ (that is in $J$ and $J_2$) for $\mathcal{H}_{RG}^{(j)}$. At second order in eccentricity, the Hamiltonian $\mathcal{H}_{RG}^{(j)}$ reads, after performing the rescaling (3.22)

$$ \mathcal{H}_{RG}^{(j)} = -\frac{9}{8}\epsilon\eta\frac{m_j}{m}\mathcal{R}_j^{-4}\left(1+\frac{4}{3}\mathcal{R}_j^{-1}X_j\bar{X}_j\right), \qquad (3.56) $$

where the quantity $\mathcal{R}_j(J)$ was defined by Eq. (3.20) and $m = \sqrt{m_1m_2}$. The total Hamiltonian that we consider for the planar relativistic co-orbital resonance is then

$$ \mathcal{H} = \mathcal{H}_K + \mathcal{H}^{(0)} + \mathcal{H}^{(2)} + \sum_{j \in \{1,2\}} \mathcal{H}_{RG}^{(j)}, \qquad (3.57) $$

where $\mathcal{H}_K$ is given by Eq. (3.24) and $\mathcal{H}^{(0)}$ and $\mathcal{H}^{(2)}$ are given by Eq. (3.29). The pair of variables $(X_j, \bar{X}_j)$ is not canonical and the equations of motion are given by Eq. (3.23). We obtain the differential system

$$ \begin{aligned} \dot{j} &= -\frac{\partial(\mathcal{H}^{(0)} + \mathcal{H}^{(2)})}{\partial \xi}, \\ \dot{\xi} &= \frac{\partial \mathcal{H}_K}{\partial J} + \frac{9}{2}\eta\epsilon \left\{ \mathcal{R}_1^{-5} \left( 1 + \frac{5}{3} \frac{X_1 \bar{X}_1}{\mathcal{R}_1} \right) - \mathcal{R}_2^{-5} \left( 1 + \frac{5}{3} \frac{X_2 \bar{X}_2}{\mathcal{R}_2} \right) \right\}, \qquad (3.58) \\ \dot{X}_j &= -2i \frac{m}{m_j} \frac{\partial \mathcal{H}^{(2)}}{\partial \bar{X}_j} + 3i\eta\epsilon \mathcal{R}_j^{-5} X_j. \end{aligned} $$

<footnote n="9">⁹A zeroth-order expansion loses all dynamical information in the circular case.</footnote>
General relativity also introduces a small perturbation on $\dot{\xi}_2$, but since the differential system does not depend on this fast angle, and we are not interested in its dynamics, we ignored the corresponding line in Eq. (3.58).

### 3.3.3 Equilibria and linearization

Taking into account the contributions of general relativity, the Lagrangian equilibria are the fixed point of the differential system (3.58). They still correspond to $X_j = 0$, and in their vicinity, the degree of freedom $(J, \xi)$ is uncoupled from the two degrees of freedom $(X_j, \bar{X}_j)$. We are once again able to study the circular and the eccentric dynamics separately.

#### Circular dynamics

We evaluate the differential system (3.58) at $X_j = 0$. Because $R_1(J = 0) = R_2(J = 0)$, and since general relativity does not perturb the value of $dJ/dt$, it is easy to see that the Lagrangian equilibria are not displaced by general relativity. The points $L_4$ and $L_5$ are still located at $(J, \xi) = (0, \pm\pi/3)$. The point $L_3$, located at $(J, \xi) = (0, \pi)$, is hyperbolic without general relativity, and given the considerations stated in Sects. 2.4.3 & 3.3.1, we know that general relativity will not make it elliptic, and we are still not interested in its dynamics.

As in Sect. 3.2.1, we linearize the system (3.58) in the vicinity of $L_4$. The resulting linear system is a small perturbation of Eq. (3.31). Writing $\Delta\xi = \xi - \pi/3$, we have

$$ \frac{d}{dt} \begin{pmatrix} \Delta \xi \\ J \end{pmatrix} = \eta \begin{pmatrix} 0 & -3 \frac{m_1 + m_2}{m} \left(1 + \frac{15}{2}\varepsilon\right) \\ \frac{9m}{4m_0} & 0 \end{pmatrix} \begin{pmatrix} \Delta \xi \\ J \end{pmatrix}. \quad (3.59) $$

The eigenvalues of this linear system can be computed immediately, as

$$ \pm i\nu \left(1 + \frac{15}{4}\varepsilon\right), \quad (3.60) $$

where $\nu$, the libration frequency of the co-orbitals without general relativity, is given by Eq. (3.32) and we recall that $v = \bar{a}\eta$. As we guessed in Sect. 3.3.1, the Lagrangian equilibria are still elliptic and general relativity is only able to slightly increase the libration frequency.

#### Eccentric dynamics

Here, we simply repeat the calculations from Sect. 3.2.2, taking into account corrections due to general relativity. At the Lagrangian point $L_4$, the eccentricities are governed by the linear differential system

$$ \frac{d}{dt} \begin{pmatrix} X_1 \\ X_2 \end{pmatrix} = \mathcal{M} \begin{pmatrix} X_1 \\ X_2 \end{pmatrix}, \quad (3.61) $$

where the matrix $\mathcal{M}$ is a perturbation of the matrix $M_0$ given in Eq. (3.41). It reads

$$ \mathcal{M} = \mathcal{M}_0 + 3i\eta\varepsilon\mathbb{I}_2 = \mathcal{M}_0 + \varepsilon\mathcal{M}_1. \quad (3.62) $$
The perturbative approach described in Sect. 2.4.2 allows the eigenvalues of $\mathcal{M}$ to be obtained without any calculation. Indeed, $\mathcal{M}_1$ is proportionnal to identity, so it keeps the same expression in any basis. Since the perturbations to the eigenvalues are given by the diagonal terms of $\mathcal{M}_1$, written in a basis that diagonalizes $\mathcal{M}_0$, we can conclude that the perturbations to the eigenvalues are simply the diagonal terms of $3i\eta\epsilon\mathbb{I}_2$. Denoting $\lambda_1$ and $\lambda_2$ the eigenvalues of $\mathcal{M}$, we thus have, at first order in $\epsilon$ (see Eq. (3.42))

$$ \lambda_1 = ig_1 + 3i\eta\epsilon, \quad \lambda_2 = ig_2 + 3i\eta\epsilon = 3i\eta\epsilon. \qquad (3.63) $$

In particular, after one orbital period, the precession of the pericentres due to general relativity, for both Lagrange and anti-Lagrange eigenmodes, amounts to

$$ \delta\omega = 6\pi\epsilon = \frac{6\pi\mathcal{G}m_0}{\bar{a}c^2}, \qquad (3.64) $$

which is the well-known expression of the pericentre advance in general relativity. The eigenvalues $\lambda_j$ are both pure imaginary, and even in the eccentric case, general relativity is unable to make the Lagrangian equilibria unstable.

## 3.4 Conclusion

In this chapter, we built an analytical model to study the planar 1 : 1 mean-motion resonance. We studied the limits of validity of the model by estimating the sizes of the remainders of the expansions that we performed. We showed that near the Lagrangian equilateral equilibria, the circular and eccentric dynamics are uncoupled, which allowed us to study them separately. By linearizing the equations of motions in the neighbourhood of the Lagrangian equilibria, we gave analytical expressions of the fundamental frequencies of the system.

In the last section of the chapter, we included the perturbations due to general relativity, showing that the relativistic effects cannot destabilize the Lagrangian equilibria. These effects are only able to slightly modify the fundamental frequencies of the system, and we analytically computed these changes.
# Co-orbital planets in a planar first-order resonance chain

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

## 4.1 Hamiltonian of the $p : p : p + 1$ resonance chain

It has been shown that, at least around low mass stars, co-orbital planets are often formed within a resonance chain (Leleu et al., 2019). In this chapter, we study an occurrence of the planar four-body problem where 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. In other words, we study the planar $p : p : p + 1$ resonance chain, where the co-orbitals have a certain nominal period and the nominal period of the outermost third planet is $(p + 1)/p$ that of the co-orbitals. We perform the calculation for any value of the integer $p$, but the figures that we show in this chapter and in Chap. 7 are for $p = 1$, that is for the resonance chain 1 : 1 : 2. We do not take into account effects due to general relativity.

### 4.1.1 The averaged Hamiltonian

We denote $(\tilde{\Lambda}_j, \tilde{D}_j; \lambda_j, -\varpi_j)$ the Poincaré canonical variables (Eq. (2.49), Sect. 2.1.4), that is, we write

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

where $\lambda_j$ and $\varpi_j$ are the mean longitude and longitude of the pericentre (of the $j^{th}$ planet), respectively. As in Chap. 3, we have defined $\beta_j = m_0 m_j / (m_0 + m_j)$ and $\mu_j = G(m_0 + m_j)$, where $G$ is the gravitational constant. According to Eq. (2.54), we write the complete Hamiltonian as

$$ H = H_K(\tilde{\Lambda}_j) + \iota H_P(\tilde{\Lambda}_j, \lambda_j, \tilde{D}_j, \varpi_j), \qquad (4.2) $$
where the Keplerian part, due to star–planet interactions, reads (Eq. (2.55))

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

and the perturbation $H_P$, due to planet–planet interaction, has a size relative to $H_K$ of order $\iota = (m_1 + m_2 + m_3) / m_0$. 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 (4.4)$$

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 $a_j^{(0)}$, and the $\tilde{\Lambda}_j$ stay close to the quantities $\Lambda_j^*$ defined as (Eq. (3.6))

$$\Lambda_j^* = m_j \sqrt{\mu_0 a_j^{(0)}}. \qquad (4.5)$$

At the exact mean motion resonance (Eq. (4.4)), the semi-major axes do not precisely verify $a_1^{(0)} = a_2^{(0)} = (p/(p+1))^{2/3} a_3^{(0)}$. However, it is still possible to consider that the exact resonance occurs at $\tilde{\Lambda}_j = \Lambda_j^*$, since Sect. 3.1.1 shows that the subsequent error is negligible with respect to the width of the resonance.

To study the dynamics in the vicinity of the resonance, the Hamiltonian is expanded in the neighbourhood of the $\Lambda_j^*$. As justified in Sect. 3.1.1, if the expansion of the Keplerian part is pushed up to the second order, then an expansion at order 0 in the perturbative part is enough, and we limit ourselves to

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

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 (4.7)$$

which is canonical if we transform the actions according to (Eq. (2.31))

$$(\tilde{\Lambda}_1, \tilde{\Lambda}_2, \tilde{\Lambda}_3, \tilde{D}_j) \mapsto (L', \Gamma', G', D'_j) = \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).$$

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. As we showed in Sect. 2.2.2, averaging over $\xi_2$ amounts in performing
a canonical transformation<fnref n="1" /> that sends the dependency on $ξ_2$ in a remainder of order $ℓ^2$, that we neglect. The averaged Hamiltonian reads

$$H'(L', \Gamma', G', \xi, D'_j, \sigma_j) = H'_K(L', \Gamma', G', D'_j) + \frac{1}{2\pi} \int_0^{2\pi} \iota H'_P(\xi, \xi_2, D'_j, \sigma_j) d\xi_2, \quad (4.8)$$

where we did not write the neglected remainder. The transformation (4.7), along with the averaging process, allows the two degrees of freedom associated with $(\Gamma', G'; \xi_2, \xi_3)$ to be eliminated, and we are left with four degrees of freedom. After the averaging process, the scaling parameter $\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 Chap. 7, it is actually more convenient to normalize by the constant (Eq. (4.5))

$$\Gamma^* = (p+1)(\Lambda_1^* + \Lambda_2^*)/p + \Lambda_3^* \quad (4.9)$$

rather than by $\Gamma'$, which is not constant when dissipation is present. In other words, we perform the canonical<fnref n="2" /> 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 (4.10)$$

while the angles are unchanged.

### 4.1.2 Expansion of the Keplerian part

As we did in Chap. 3, we expand the Keplerian part of the Hamiltonian 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 (4.11)$$

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 (4.12)$$

The third term is constant and can be removed without changing the dynamics. In the variables $L, G, \Gamma$ and $D_j$, $\mathcal{H}_K$ takes the form

$$\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 (4.13)$$

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 (4.14)$$

<footnote n="1">¹For convenience, we do not change the name of the variables.</footnote>

<footnote n="2">²The Jacobian of this transformation is not symplectic, but the rescaling of the Hamiltonian guarantees that the form of Eq. (2.19) is preserved.</footnote>
Without dissipation, $\Gamma$ and $G$ are both integrals of motion. Due to the normalization by the constant $\Gamma^*$, $\Gamma$ is close to 1 and we simply evaluate $\mathcal{H}_K$ at $\Gamma = 1$, hence achieving the reduction to only one parameter. The relevant parameter to consider is $\Gamma/G$ and the variations of $\Gamma$ are now<fnref n="3" /> absorbed by $G$ (Eqs. (4.34) and (4.35)). The last term of Eq. (4.13) is constant and can be removed from the Hamiltonian.

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 (4.5). In this case the approximation $\Gamma = 1$ becomes $\Delta \Gamma = 0$ and Eq. (4.11) yields, once normalized,

$$ \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 (4.15) $$

where $\Delta\Upsilon = \Delta G + D_1 + D_2 + D_3 = \sum_j \Delta\Lambda_j$. We find the Hamiltonian (4.15) to be well adapted to the analytical work derived in Sect. 4.3.1, while we use the Hamiltonian (4.13) elsewhere. Moreover, we do not perform the evaluation $\Gamma = 1$ in Chap. 7, where we study the resonance chain $p : p : p+1$ in presence of tidal dissipation, making $\Gamma$ a variable quantity. Hamiltonians (4.13) & (4.15) both yield the same dynamics and using either is only a matter of preference.

### 4.1.3 Expansion of the perturbative part

The perturbative part $\iota \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

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

We denote $X_j = \sqrt{2\tilde{D}_j/\tilde{\Lambda}_j} e^{i\tilde{\omega}_j} = e_j e^{i\tilde{\omega}} + \mathcal{O}(e_j^3)$. For a pair $(p_1, p_2) \in \{(1, 2), (1, 3), (2, 3)\}$ of planets, the Hamiltonian $\mathcal{H}_{p_1,p_2}$ is given by Eq. (2.58). The d'Alembert rule (2.59), combined with the average over $\xi_2 = p\lambda_2 - p\lambda_3$, implies that $\mathcal{H}_{1,2}$ has no term of odd degree in eccentricity, while $\mathcal{H}_{1,3}$ and $\mathcal{H}_{2,3}$ have no term of degree 0. Since we limit ourselves to the second degree 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 (4.17) $$

where $j \in \{1, 2\}$. The superscript refers to the degree in eccentricity, while the subscript refers to the considered pair of planets. The Hamiltonian $\mathcal{H}^{(0)}$ has no subscript since only the co-orbital pair yields terms of degree 0 and no confusion is possible. The part $\mathcal{H}_{1,2}$ of the Hamiltonian was already computed in Chap. 3, and is given by Eq. (3.29). Using the variables $D_j$ and $\sigma_j$, and given the normalization (4.10), we can write

$$ \begin{aligned} \mathcal{H}^{(0)} &= \frac{m_1}{m_0 C_2} (\cos \xi - \Delta^{-1}), \\ \mathcal{H}_{1,2}^{(2)} &= \frac{m_1}{m_0 C_2} \left\{ A_h (C_1 D_1 + C_2 D_2) + 2\sqrt{C_1 C_2 D_1 D_2} \operatorname{Re} (B_h e^{i(\sigma_2 - \sigma_1)}) \right\}, \end{aligned} \quad (4.18) $$

where $\Delta = \sqrt{2 - 2 \cos \xi}$. The coefficients $A_h$ and $B_h$ are given by Eq. (3.30), and the $C_j$ that depends only on the planetary masses and on $p$, were defined by Eq. (4.14).

<footnote n="3">³After the evaluation $\Gamma = 1$, $G$ is in fact $G/\Gamma$, but we abusively keep writing $G$.</footnote>
Following Laskar and Robutel (1995) and using the algebraic manipulator Trip (Gastineau and Laskar, 2011), $\mathcal{H}_{j,3}^{(1)}$ and $\mathcal{H}_{j,3}^{(2)}$ 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(p\delta_{j,1}\xi - \sigma_j) + C_{p,2}^{(1)} \sqrt{2C_3 D_3} \cos(p\delta_{j,1}\xi - \sigma_3) \right\}, \quad (4.19) $$

and

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

where $\delta_{j,1} = 1$ if $j = 1$ and zero if $j = 2$. 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 Eqs. (B.2) & (B.3) of Appendix B. The final Hamiltonian of the model is then

$$ \mathcal{H}(L, D_j; \xi, \sigma_j) = \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 (4.21) $$

As opposition to the complete Hamiltonian of Eq. (4.2), we refer to Eq. (4.21) as the simplified Hamiltonian, or as the model. We recall the three simplifications that led from Eq. (4.2) to Eq. (4.21)

* Average over the angle $\lambda_2 - \lambda_3$, that is, Lie serie expansion at first order in $\iota$, while ignoring the associated near-identity change of variables (see Sect. 2.2.2).

* Expansion at second order in eccentricity.

* Expansion of the semimajor axes in the neighbourhood of Eq. (4.5), at second order for $H_K$, and at the zeroth order for $\iota H_P$.

We denote $F_0 : \mathbb{R}^8 \to \mathbb{R}^8$ the vector field derived from Eq. (4.21) by the Hamilton Eqs. (2.18).

# 4.2 Important theoretical results on resonance chains

## 4.2.1 From fixed points to libration centres

As we will see in Sect. 4.3.2, the model (4.21) has equilibria, that is, points in the phase space where the vector field $F_0$ vanishes. 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 exists constants $c_j$ such that

$$ \sigma_j = -p\lambda_2 + (p+1)\lambda_3 - \varpi_j = c_j. \quad (4.22) $$

However, the secular angle $\xi_3 = -p\lambda_2 + (p+1)\lambda_3$ and the pericentres $\varpi_j$ are not constant at the equilibria, but they all precess with the same frequency, which we denote $\nu_3$.
Even though the model has equilibria, the complete Hamiltonian (4.2) has none, since no initial condition in the n-body problem leads to the immobility of the system. The canonical change of variable that we ignored, inherent to the Lie serie expansion, explains the inconsistency. Indeed, as is explained in Sect. 2.2.2, the average over the angle $\xi_2$ in Eq. (4.8) is in fact a near-identity canonical transformation that allows the dependency on that angle to be pushed to the neglected second order in $\iota$. As a consequence, the simplified Hamiltonian $\mathcal{H}$ (Eq. (4.21)) does not depend on $\boldsymbol{x} = {}^t(L, \Gamma, G, D_j; \xi, \xi_2, \xi_3, \sigma_j)$, but rather on $\tilde{\boldsymbol{x}} = {}^t(\tilde{L}, \tilde{\Gamma}, \tilde{G}, \tilde{D}_j; \tilde{\xi}, \tilde{\xi}_2, \tilde{\xi}_3, \tilde{\sigma}_j)$, where (Eq. (2.63))

$$ \tilde{\boldsymbol{x}} = \Phi_{\iota\chi}(-1, \boldsymbol{x}) = e^{-\iota L_{\chi}} \boldsymbol{x}. \qquad (4.23) $$

Equations. (2.67) to (2.70) show that the transformation $\Phi_{\iota\chi}(-1, \cdot)$ is quasiperiodic with vector of frequencies $\omega_0 = \partial\mathcal{H}_K/\partial\boldsymbol{J}$, where $\boldsymbol{J} = {}^t(L, \Gamma, G, D_1, D_2, D_3)$. Furthermore, at an equilibrium of the model, because $\tilde{\boldsymbol{x}}$ is constant, the Fourier decomposition of $\chi$ in Eq. (2.69) only depends on the 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}, \qquad (4.24) $$

where these quantities are evaluated at the equilibrium. This means that at an equilibrium of the model (4.21), the transformation $\tilde{\boldsymbol{x}} \mapsto \boldsymbol{x}$ is quasiperiodic with frequencies $\nu_2$ and $\nu_3$, and so is $\boldsymbol{x}$. More precisely, $\boldsymbol{x}$ is periodic with frequency $\nu_2$, in a frame rotating with the pericentres $\varpi_j$ at frequency $\nu_3$. Furthermore, for a quantity invariant around the axis of total angular momentum, like $\xi$ or the $\sigma_j$, the Fourier decomposition does not feature $\nu_3$, and the motion of $\boldsymbol{x}$ in the complete Hamiltonian (4.2) is periodic with frequency $\nu_2$.

In other words, the motion of a quantity $z$ depending on $\boldsymbol{x}$, along the flow of the complete Hamiltonian (Eq. (4.2)), using an equilibrium of the simplified Hamiltonian (Eq. (4.21)) as initial condition, is

* Periodic with frequency $\nu_2$ if $z$ is invariant by the angle translation $(\lambda_j, \varpi_j) \mapsto (\lambda_j + \vartheta, \varpi_j + \vartheta)$, where $\vartheta$ is any constant angle.

* Quasiperiodic with frequencies $\nu_2$ and $\nu_3$ if $z$ is not invariant by this angle translation.

Given the general change of variables proposed by Eq. (A.1) of Delisle (2017), it is easy to notice that this result is true for any chain of resonances. 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.

### 4.2.2 An algorithm to find the libration centres

According to Sect. 4.2.1, the equilibria of the simplified Hamiltonian (Eq. (4.21)), should correspond to points in the phase space where, along the flow of the complete Hamiltonian (Eq. (4.2)), the motion of a quantity $z$ invariant by rotation and depending on $\boldsymbol{x} = {}^t(L, \Gamma, G, D_j; \xi, \xi_2, \xi_3, \sigma_j)$, is periodic with frequency $\nu_2$. However, due to the expansion of the Hamiltonian in eccentricities and semimajor axes, the *libration centres* of the complete Hamiltonian do not coincide with the equilibria of the model. We develop here an iterative algorithm, similar to that in Couetdic et al. (2010), to find a *libration*
<figure><img src="image_7.png" /><figcaption>Fig. 4.1 — Trajectories of $e_j e^{i\sigma_j}$ in the complete Hamiltonian (4.2) for the six iterations needed for the algorithm to converge to the *libration centre*. Iteration 0 is an elliptic equilibrium$^4$ of the simplified Hamiltonian (4.21). It is rather far from the *libration centre* of the complete Hamiltonian (see also Fig. 4.4). After six iterations, the algorithm converges to the *libration centre* and the motion is periodic (hence the closed curves) with frequency $\nu_2$ (Sect. 4.2.1). The planetary masses are those of Fig. 4.2 and the resonance chain is 1 : 1 : 2.</figcaption></figure>

*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 (4.2), that the trajectories are quasiperiodic, 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 (4.25)$$

where the coordinates of the vector $\omega = t(\nu, \nu_2, \nu_3, g_1, g_2, g_3)$ are the fundamental frequencies of the complete Hamiltonian (4.2). The frequencies $\nu, \nu_2$, and $\nu_3$ have approximate values given by Eqs. (4.39), (4.24), and (4.30), respectively. The description of the algorithm is as follows:

* Find the position of an elliptic equilibrium of the simplified Hamiltonian (4.21) with a Newton-Raphson method. Use it as the initial condition to integrate numerically the trajectories of the complete Hamiltonian (4.2).

* For a complex quantity $z$ depending on the trajectories of (4.2), obtain the decomposition (4.25) using a frequency analysis method (e.g. Laskar, 1993).

<footnote>$^4$It is the main branch, introduced in Sect. 4.3.1, at $\delta = 5$ (see Eq. (4.34)).</footnote>
• Identify the terms depending on frequencies other than $v_2$ and set to 0 the corresponding coefficient $z_k$.

• Proceed similarly for enough 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 (4.21).

The process is iterated until a convergence occurs. In Fig. 4.1, 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 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$, instead of setting to 0 the coefficients $z_k$ of the terms that depends on frequencies other than $v_2$. This does not prevent the algorithm from properly converging.

## 4.3 Equilibria and linearization

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

### 4.3.1 Analytical results

The fixed points of the Hamiltonian (4.21) cannot be given analytically, even if it is truncated at first degree in eccentricity. Similar difficulties were met by Delisle (2017) for resonance chains with first-order resonances between non-consecutive planets. However, we show here that further simplifications allow us to obtain analytical expressions of the equilibria and of the eigenvalues of the linearized system. In this subsection, we consider the Hamiltonian $\mathcal{H} = \mathcal{H}_K + \mathcal{H}^{(0)} + \mathcal{H}^{(1)}$, that is, the Hamiltonian of Eq. (4.21), where the terms of second order in eccentricity (those featuring the superscript (2)) have been removed. The Hamiltonian $\mathcal{H}_K + \mathcal{H}^{(0)} + \mathcal{H}^{(1)}$ is only a small perturbation of the Hamiltonian $\mathcal{H}_K + \mathcal{H}^{(0)}$, studied in Sect. 3.2.1, which has equilibria at $\xi \in \{\pm\pi/3, \pi\}$. The equilibrium at $\xi = \pi$ is unstable and we are not interested in its dynamics. Only the equilibria $\xi = \pm\pi/3$ are interesting to us, and since their dynamics are symmetric, it is enough to only consider $\xi = \pi/3$.

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$). We proceed in two steps :

• We first evaluate $\mathcal{H}^{(1)} = \mathcal{H}_{1,3}^{(1)} + \mathcal{H}_{2,3}^{(1)}$ at $\xi = \pi/3$, meaning that $\mathcal{H}^{(1)}$ loses its dependency on $\xi$.

• We then replace the variable $\Delta L$ by the constant $\Delta L^*$ in the anti-diagonal term<fnref n="5" /> of $\mathcal{H}_K$ in (4.15), where $\Delta L^*$ is the value of $\Delta L$ for which $\partial\mathcal{H}_K/\partial\Delta L$ vanishes, that is

<footnote n="5">⁵The term proportional to $\Delta L \Delta \Upsilon$</footnote>
$$ \Delta L^* = -\frac{pC_2}{C_1 + C_2} (\Delta G + D_{1,0} + D_{2,0} + D_{3,0}). \quad (4.26) $$

The $D_{j,0}$ are given by Eq. (4.29). While the evaluation at $\xi = \pi/3$ allows analytical expressions for the position of the fixed points, the evaluation at $\Delta L = \Delta L^*$ uncouples $(\xi, \Delta L)$ from $(\sigma_j, D_j)$ and enables analytical expressions of the eigenvalues of the linearized system in the vicinity of the fixed points. The differential system derived from $\mathcal{H}_K + \mathcal{H}^{(0)} + \mathcal{H}^{(1)}$, after 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 \sigma_{3,0} = \arctan \frac{m_1 \sin p\xi_0}{m_2 + m_1 \cos p\xi_0} + (1 - \epsilon)\pi, \quad (4.27) $$

where<fnref n="6" />

$$ \epsilon = \begin{cases} 0 & \text{if } \nu_3 < 0, \\ 1 & \text{if } \nu_3 > 0, \end{cases} \quad (4.28) $$

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

$$ \begin{gathered} C_1 D_{1,0} = C_2 D_{2,0}, \quad \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 \quad \text{and} \\ \frac{C_1 C_{p,1}^{(1)2} m_1^2 p^2}{2 C_3^2 m_0^2 (p+1)^2} = \left(\frac{\nu_3}{\eta}\right)^2 D_{1,0} \end{gathered} \quad (4.29) $$

where the precession frequency of the pericentres is

$$ \nu_3 = -\eta K (\Delta G + D_{1,0} + D_{2,0} + D_{3,0}), \quad K = \frac{3p^2 C_1 C_2}{C_1 + C_2} + 3p(p+1)C_3, \quad (4.30) $$

and we define the constant $H$ by

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

The unknowns quantities of Eq. (4.29) 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 quantity $D_{1,0}$. The precession frequency of the pericentres, $\nu_3$, depends on $D_{1,0}$ (see Eq. (4.30)). 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$, we can rewrite Eq. (4.29) as a third-degree polynomial in $Z$

$$ \begin{gathered} Z^3 - PZ - Q = 0, \\ \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^2 C_3^2 K^2 m_0^2 (p+1)^2}. \end{gathered} \quad (4.32) $$

The coefficients $P$ and $Q$ of this polynomial depend on the parameter $\Delta G$, that we properly define in Eq. (4.35). There is a bifurcation between 1 and 3 real solutions when $27Q^2 - 4P^3 = 0$ (Cardano, 1545), 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}. \quad (4.33) $$

<footnote n="6">⁶An analytical expression for $\nu_3$ is given by Eq. (4.30).</footnote>

<footnote n="7">⁷$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>
<figure><img src="image_8.png" /><figcaption>Fig. 4.2 — 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. (4.32). Right: Values of $\nu/\eta$ and $\nu_3/\eta$ along the main branch as a function of $\delta$, for the same resonance chain, predicted by Eqs. (4.39) and (4.30). 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. (4.29), $e_1 = e_2$, and for this choice of masses and resonance chain $e_3/e_1 = 0.5471$. In Table 4.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 $\nu$ and $\nu_3$ are shown on the right, and are also easy to spot in Fig. 4.6.</figcaption></figure>

In the rest of this work we use the parameter

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

where

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

In this chapter, $\Gamma = \Gamma'/\Gamma^* \approx 1$ is simply evaluated at 1 and ignored, but not in Chap. 7, 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 and Lemaitre (1983). The fixed points are given by the roots of the third-degree polynomial in $Z$ (4.32), which has one or three real solutions depending on $\delta$, hence a bifurcation. The solutions of Eq. (4.32) are plotted in Fig. 4.2. 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. 4.3.2 that the topology of the Hamiltonian (4.21) is different (see Table 4.1). However, we show in Fig. 4.3 that these analytical expressions are accurate for low eccentricities.

We now linearize the set of Eqs. (C.1) in the vicinity of the main branch. We make
use of the cartesian coordinates

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

and denoting $X = {}^t(u_1, u_2, u_3, v_1, v_2, v_3, \Delta L, \xi)$, the linearized system is

$$\begin{gathered} \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{gathered} \qquad (4.37)$$

and $X_0$ is the equilibrium value of $X$. The matrix $Q_6$ is given by Eq. C.4. Its characteristic polynomial reads

$$\det(\lambda \mathbb{I}_6 - Q_6) = (\lambda^2 + \nu_3^2)^2 \left( \lambda^2 + \nu_3^2 - 2\nu_3 I \sum_{j=1}^{3} D_{j,0} \right), \qquad (4.38)$$

where the constant $I$ is defined in Eq. (C.3). It is interesting to note that the precession frequency of the pericentres, $\pm i\nu_3$, is an eigenvalue of the differential system (4.37). 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 (see their Eq. (109)). The eigenvalues of $Q_2$ are $\pm i\nu$, where $\nu$ is given by Eq. (3.32)

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

We recall that $\nu$ is the libration frequency of the angle $\xi$ in the neighbourhood of the equilateral Lagrangian configuration. Figures 4.2 and 4.5 show that $\nu_3 < 0$ for the main branch, which ensures that the roots of (4.38) are pure imaginary. Evaluating the eigenvalues $\pm i\nu$ and $\pm i\nu_3$ in the vicinity of the main branch shows that, at $\delta \approx -5.6$ for the planetary masses in Fig. 4.2, 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 $\nu$ and $\nu_3$ are shown in Fig. 4.2, and are also very visible on the stability map from Fig. 4.6. We show in Sect. 7.2.2 that the secular resonance 1 : 1 has important consequences for the tidal stability of the co-orbital pair.

## 4.3.2 Numerical results

In the previous section, strong approximations allowed the equilibria and eigenvalues of the model to be given analytically. We now develop a Newton-Raphson-based algorithm to numerically and exactly find the equilibria of the Hamiltonian (4.21) and compare them with the analytical approximations elaborated in Sect. 4.3.1.

The Newton-Raphson algorithm finds equilibria by starting from a random initial condition in the phase space, and we discretize it by exploring a grid in the parallelepiped defined by $|u_j| < 0.08$ and $|v_j| < 0.08$ (recall that $e_j^2 = C_j(u_j^2 + v_j^2)$). In order to pretend to have found all possible equilibria, it is necessary that the discretization is thin, or in other
<table><caption>Table 4.1 — The 15 equilibria of the simplified Hamiltonian (4.21). They 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. 4.2. Branch 2 is hyperbolic only for $5.548 \le δ \le 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}}$ (It occurs at $δ = 1$ in the analytical results of Sect. 4.3.1). Branch 1 is the only elliptic branch existing for all values of $δ$ and it is the main branch introduced in Sect. 4.3.1. It corresponds to the only real solution of Eq. (4.32) 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 (4.2) has more libration centres, either because we did not discretize the phase space thinly enough to find them or because they do not exist at second order in eccentricity.</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>

words, that the grid is dense. However, since the phase space has eight dimensions, a thin discretization leads to too many initial conditions, and to limit the computational time, we are only able to discretize the axes $u_j$ and $v_j$ with eight points, evenly distributed in the range $[-0.08, 0.08]$. Since the semimajor axes stay close to their Keplerian value, the L-axis is discretized with the unique point $L_0 = L^* := \Lambda_1^*/\Gamma^*$, which corresponds to the value of L at the Keplerian resonance (Eq. (4.5)). Similarly, because the model (4.21) is only a perturbation of $H_K + H^{(0)}$, its equilibria are expected to be close to $ξ_0 \in \{\pm\pi/3, \pi\}$ (see Fig. 3.1) and the ξ-axis is discretized with these three points only. The algorithm hence explores a grid composed of $3 \times 8^6$ points in the phase space.

For $δ = 7$, all 786 432 initial conditions of the Newton–Raphson method converged towards 15 distinct equilibria. The vector field $F_0$ derived from the Hamiltonian (4.21) 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. In Table 4.1, we display the 15 found equilibria, their hyperbolic or elliptic nature, and the value of $δ$ that gives birth to the branch. Since the Hamiltonian (4.21) 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 4.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 4.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. 4.3 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, we also plot them as given by Eqs. (4.27) and (4.32), for comparison. This section shows how the analytical results are unable
<figure><img src="image_9.png" /></figure>

Fig. 4.3 — Position of the elliptic branches 1, 2, and 3 (see Table 4.1) of equilibria of the Hamiltonian (4.21) 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 Fig. 4.2. 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. (4.27) and (4.32), is plotted by thin grey lines. 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. (4.27). 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.

to locate the equilibria of the model (4.21) for $e_j \gtrsim 0.005$ (see Fig. 4.3) 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 4.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.

# 4.4 Comparison with the unaveraged system

## 4.4.1 Fixed points versus libration centres

As previously justified, the *libration centres* of the complete Hamiltonian (4.2) (Sect. 4.2.1) do not coincide with the equilibria of the model (4.21). We use the algorithm described in Sect. 4.2.2 to locate the *libration centres* of the complete system and compare them with the equilibria of the model. Since this algorithm is only able to find *libration centres* associated with elliptic fixed points, we locate the branches of *libration centres* associated with branches 1, 2, and 3 of Table 4.1. We confirm the existence of a small hyperbolic zone for branch 2 in the complete Hamiltonian when the algorithm stops converging as it
<figure><img src="image_10.png" /><figcaption>Fig. 4.4 — Position of branch 1 of elliptic libration centres of the complete Hamiltonian (4.2) in the resonance 1 : 1 : 2 for $-7 \le \delta \le 7$. The planetary masses are the same as in Fig. 4.2. As a comparison, branch 1 of equilibria of the simplified Hamiltonian (4.21) 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. Then, the simplified and the complete Hamiltonian diverge when $\delta \to +\infty$.</figcaption></figure>

travels along the branch (by slowly incrementing the value of $\delta$). We plot in Fig. 4.4 the main branch (branch 1) of the libration centres of the complete Hamiltonian (4.2) and we compare it with the main branch of equilibria of the simplified Hamiltonian (4.21). In order not to overload the figure, only branch 1 is plotted. It is also the most interesting branch when tidal dissipation is involved (see section 7.2).

In Fig. 4.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 (4.30). 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.

### 4.4.2 A 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. 4.2) 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
<figure><img src="image_11.png" /><figcaption>Fig. 4.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.4. The analytical expression (4.30) is plotted in grey for comparison. The colour gives the value of $n_1/n_3$. As expected from Eqs. (4.30) & (4.34), $v_3$ is proportional to $\delta$ far from the resonance, where the eccentricities are small.</figcaption></figure>

of the algorithm described in Sect. 4.2.2, 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 and 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 $\nu$ is extracted from the frequency analysis of $e^{i\xi}$. We obtain two values of $\nu$, namely $\nu^{(1)}$ for the first half of 40 000 periods and $\nu^{(2)}$ for the second half. The diffusion index, defined as (Robutel and Gabern, 2006)

$$ \zeta_{\nu} = \log_{10} \left| \frac{\nu^{(1)} - \nu^{(2)}}{\nu^{(1)}} \right|, \qquad (4.40) $$

measures the degree of quasi-periodicity of the orbit. Orbits with $\zeta_{\nu} < -6$ are considered close to be quasi-periodic (stable), while orbits with $\zeta_{\nu} > -2$ are very chaotic (unstable). We plot the stability map for the resonance chain 1 : 1 : 2 in Fig. 4.6. Secular resonances between $\nu$ and $v_3$, already predicted by the analytical results in Fig. 4.2 are visible and induce chaotic motion. Overall, this stability map shows that the chain $p : p : p + 1$ is mainly stable.

## 4.5 Conclusion

In this chapter, we have studied the dynamics of a pair of co-orbital planets in the presence of a first-order mean motion resonance with a third planet orbiting outside the co-orbitals.
<figure><img src="image_12.png" /><figcaption>Figure 4.6 — Diffusion index $\zeta_{\nu}$ as a function of $\delta$ and $\xi_{0}$. The planetary masses are as in Fig. 4.2 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 $\nu$ and the precession frequency of the pericentres $\nu_{3}$, predicted by the analytical results (see Fig. 4.2), 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 $\nu$ and $\nu_{3}$ is tidally attractive (see Fig. 7.1). Systems undergoing tides entering this zone can either converge towards the top of the figure or become completely chaotic.</figcaption></figure>

We built an averaged model (Eq. (4.21)) from which we were able to analytically retrieve the equilibria of the system and the eigenvalues in their vicinity. With these analytical results, we showed that the resonance chain $p : p : p + 1$ has some similarities with the two-planet second fundamental model of resonance of Henrard and Lemaitre (1983), as well as some other proper characteristics, like the existence of a secular 1 : 1 resonance between the libration frequency of the co-orbitals and the precession frequency of the pericentres, that can lead to chaotic orbits.

The analytical expressions came from strong assumptions, and we have compared them with numerical results, obtained after all simplifications were removed. This allowed us to reach the conclusion that the analytical expressions were correct, for eccentricities less than 0.005. We have shown that the resonance chain $p : p : p + 1$ possesses many families of equilibria, some of them being stable. We also formally studied the relation between the equilibria of the averaged model and the *libration centres* of the unaveraged, complete problem, and we have checked the ability of our model to predict the location of the *libration centres* of the system, by means of a dedicated algorithm.
# Main theory on tidal dissipation

This chapter aims to establish the main theory of tides, that we use in Chaps. 6 and 7.

## 5.1 Tidal potentials

Celestial bodies are not perfectly rigid and the differential gravitational interactions between them give rise to inelastic deformations and distortions, leading to dissipation of energy by heating up. This energy loss modifies their spins and orbits and it is important to have a satisfying tidal description to study the long term dynamics of planetary systems. Our current mathematical understanding of tidal mechanisms comes from a very general formulation initiated by Darwin (1880). We consider in this chapter an extended body $\mathbb{B}$ of barycentre $O$ and mean radius $R_0$ upon which tides are raised.

### 5.1.1 The perturbing potential

Tides arise when a celestial body, referred to as perturbing body, orbits in space around $\mathbb{B}$. We note $\tilde{\mathbf{r}}$ the position of the perturbing body with respect to $O$. If $\mathbf{r}$ is the position, with respect to $O$, of a point within $\mathbb{B}$, this point undergoes from the perturbing body the potential per unit mass

$$W(\mathbf{r}) = -\frac{G\tilde{m}}{|\mathbf{r} - \tilde{\mathbf{r}}|} = -\frac{G\tilde{m}}{\tilde{r}} \sum_{l=0}^{+\infty} \left(\frac{r}{\tilde{r}}\right)^l P_l(\cos S), \quad (5.1)$$

where $\tilde{m}$ is the mass of the perturbing body, $S$ is the angle between $\mathbf{r}$ and $\tilde{\mathbf{r}}$ and the $P_l$ are the Legendre polynomials, whose definition is recalled in appendix A. Tides are due to the differential acceleration caused by $W$ within $\mathbb{B}$, which means that the first two terms of the perturbing potential, those corresponding to $l = 0$ and $l = 1$, do not contribute to the deformation of $\mathbb{B}$. Indeed, their gradient complies with

$$\nabla \left( \frac{G\tilde{m}}{\tilde{r}} \left( 1 + \frac{\mathbf{r} \cdot \tilde{\mathbf{r}}}{\tilde{r}^2} \right) \right) = \frac{G\tilde{m}\tilde{\mathbf{r}}}{\tilde{r}^3}, \quad (5.2)$$
and these terms are merely responsible for the orbital motion of $\mathbb{B}$ around the perturbing body. Hence we redefine the perturbing potential as

$$W(\mathbf{r}) = -\frac{\mathcal{G}\tilde{m}}{\tilde{r}} \sum_{l=2}^{+\infty} \left(\frac{r}{\tilde{r}}\right)^l P_l(\cos S). \qquad (5.3)$$

In most practical cases, $r \ll \tilde{r}$, and only the term $l=2$ is retained

$$W(\mathbf{r}) = -\frac{\mathcal{G}\tilde{m}r^2}{\tilde{r}^3} P_2(\cos S). \qquad (5.4)$$

### 5.1.2 Perturbed potential and definition of the Love numbers

In order to study the evolution of a planetary system undergoing tidal dissipation, we need to establish a relation between the perturbing potential $W(\mathbf{r}, t)$ and the perturbed potential $V(\mathbf{r}, t)$, that is, the potential raised by the redistribution of mass inside $\mathbb{B}$. Generally, the auto-gravitation is the main contribution to the equilibrium figure of a body, and tides raised by other bodies of the system act as tiny perturbations. Thus, it is safe to assume that $V$ responds linearly to $W$ and we write $V = K(W)$, where $K$ is a linear operator. Furthermore, $V(\mathbf{r}, t)$ should a priori depend on the value of $W$ at any point in space and at any time prior to $t$. To account for this and for the linearity, we express $V(\mathbf{r}, t)$ as the convolution product (Boué et al., 2019)

$$V(\mathbf{r}, t) = \int_{-\infty}^{t} \int_{\mathbb{R}^3} k(\mathbf{r}, \mathbf{r}', t-t') W(\mathbf{r}', t') d\mathbf{r}' dt'. \qquad (5.5)$$

The integration over space can be restricted to the volume of $\mathbb{B}$, and in fact, when $\mathbb{B}$ is incompressible, it can even be restricted to any sphere encircling completely the body $\mathbb{B}^1$. For simplicity, we choose a sphere $\mathcal{S}$, centered on $\mathcal{O}$, and of radius $R$ slightly larger than $R_0$, as to ensure that $\mathbb{B}$ is entirely inside of $\mathcal{S}$. Thus, we write

$$V(R, \sigma, t) = \int_{-\infty}^{t} \int_{\mathcal{S}} k(R, \sigma, \sigma', t-t') W(R, \sigma', t') d\sigma' dt', \qquad (5.6)$$

where $\sigma = (\theta, \varphi)$ are the colatitude and longitude and $d\sigma = \sin\theta d\theta d\varphi$. We now decompose $V$ and $W$ over the spherical harmonics (Appendix A.1), that is

$$\begin{aligned} V(R, \sigma, t) &= \sum_{l=0}^{+\infty} \sum_{m=-l}^{l} V_{lm}(R, t) Y_{lm}(\sigma), \\ W(R, \sigma', t') &= \sum_{l'=2}^{+\infty} \sum_{m'=-l'}^{l'} W_{l'm'}(R, t') Y_{l'm'}(\sigma'), \end{aligned} \qquad (5.7)$$

where $l'$ starts at 2, as in Eq. (5.3). When $\mathbb{B}$ is isotropic, the coefficients of the decomposition verify a relation similar to Eq. (5.6), namely

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

<footnote n="1">¹This comes from the general solution of the Laplace equation in spherical coordinates (Eq. (A.11)), that forces the form of $W$ everywhere inside the sphere, once $W$ is known on the surface of the sphere.</footnote>
The demonstration of Eq. (5.8) (Boué et al., 2019) is given in appendix A.3, while we recall the definition of $Y_{lm}$ in appendix A.1. In this equation, the Love numbers $k_{lm}$ depend on the arbitrary radius $R$ and we would like to define more general Love numbers. According to the Poisson equation, we have $\Delta V = 0$ everywhere outside of $\mathbb{B}$ and $\Delta W = 0$ everywhere in space (expect inside the body responsible for the tides). Since $V$ and $W$ do not diverge at $\mathcal{O}$ and at infinity, we can deduce, from the general solution of the Laplace equation $\Delta = 0$ in spherical coordinates (Eq. (A.11)), that

$$V_{lm}(R, t) \propto R^{-(l+1)} \quad \text{and} \quad W_{lm}(R, t) \propto R^l, \qquad (5.9)$$

in the vicinity of $\mathbb{S}$. As a consequence, we have

$$k_{lm}(R, t) \propto R^{-(2l+1)}, \qquad (5.10)$$

and we redefine the Love numbers as

$$k_{lm}(t) = \left(\frac{R}{R_0}\right)^{2l+1} k_{lm}(R, t). \qquad (5.11)$$

We also redefine the coefficients $V_{lm}$ with an equation similar to (5.8)

$$V_{lm}(t) = \int_{-\infty}^{t} k_{lm}(t-t')W_{lm}(R_0, t')dt', \qquad (5.12)$$

which yields, for $r \ge R$

$$V(r, \sigma, t) = \sum_{l=2}^{+\infty} \left(\frac{R_0}{r}\right)^{l+1} \sum_{m=-l}^{l} V_{lm}(t)Y_{lm}(\sigma). \qquad (5.13)$$

Equations (5.12) and (5.13) give the final form of the perturbed potential, as a function of the perturbing potential. To give an expression for $V$, the Love numbers $k_{lm}(t)$ need to be known. Choosing a particular tidal model is equivalent to choosing a particular expression for the $k_{lm}(t)$. The convolution product (5.12) can be conveniently written in the frequency domain

$$\hat{V}_{lm}(\eta) = \hat{k}_{lm}(\eta)\hat{W}_{lm}(R_0, \eta), \qquad (5.14)$$

where $\hat{f}$ denotes the Fourier transform of $f$. In the rest of this manuscript, we drop the hat to denote the Fourier transform. In general, $R_0/r$ is much smaller than unity and the serie (5.13) is truncated at $l = 2$, which means that the knowlegde of $k_{2m}(t)$, for $-2 \le m \le 2$, is enough. Furthermore, the dependency on $m$ is artificial, since there exists integers $p$ and $q$ such that $k_{2m}(\eta) = k_{20}(pn + q\omega)$, where $\omega$ and $n$ are the rotation rate and mean motion of the perturbed body (around the perturbing body), respectively (e.g. Correia and Laskar, 2003, Sect. 2.1). Only the second Love number $k_2(t) := k_{20}(t)$ is thus needed to define a tidal model.

## 5.2 Simple rheologies and subsequent tidal models

In order to define a tidal model, we need to obtain an expression for the Love number $k_2(t)$. The idea is to consider a case where the perturbing potential $W$ and the perturbed
potential $V$ are both independently known and to use Eq. (5.14) to deduce the $k_{lm}$. Since these numbers depend on the considered body, but not on the choice of potential, this method allows us to determine them in a convenient way.

We consider for $W$ the potential generated by the own rotation of $\mathbb{B}$. This rotation is responsible for an equatorial bulge whose associated potential is $V$. These potentials are axisymmetric and the coefficients of their decomposition (5.7) will verify $W_{lm} = V_{lm} = 0$ for $m \neq 0$. Hence, the decomposition in the spherical harmonics will conveniently feature only the Legendre polynomials $P_l$, since $Y_{l0} \propto P_l$ (appendix A). Besides being axisymmetric, $V$ and $W$ are also independent of time, which means that $V$ depends on $W$ only at present time. Given the form of Eq. (5.12), this implies $k_2(t) \propto \delta(t)$ where $\delta$ is Dirac's distribution. This choice for the potentials implies a frequency of excitation zero, and we will only retrieve the value of the Love number $k_2$ at frequency $\eta = 0$. Indeed, $k_2(t) \propto \delta(t)$ yields $k_2(\eta) = C_{\text{stant}} = k_2(\eta = 0)$. To obtain the true dependency of $k_2$ on the frequency, potentials more general than $V$ and $W$ need to be chosen. In Sect. 5.2.2, we present the general method to retrieve the dependency on $\eta$.

### 5.2.1 The elastic body

In this section, $\mathbb{B}$ is homogeneous, incompressible and perfectly elastic with a shear modulus $\mu$. We note $m$ and $\omega$ the mass and rotation rate of $\mathbb{B}$ and define $\omega_c$ as the rotation rate that would compensate gravity at the equator if $\mathbb{B}$ was rigid, that is, $\omega_c^2 = Gm/R_0^3$. We also introduce the small parameter $\iota = \omega^2/\omega_c^2$. With these notations, the centrifugal acceleration in the frame rotating with $\mathbb{B}$ derives from the potential

$$W(r, \theta) = -\frac{1}{2}\iota\omega_c^2 r^2 \sin^2 \theta = \frac{1}{3}\iota\omega_c^2 r^2 [P_2(\cos \theta) - 1]. \quad (5.15)$$

To determine the potential $V$, We consider that the surface of $\mathbb{B}$ is deformed axisymmetrically under its own rotation (e.g. Ragazzo and Ruiz, 2015; Wahr, 1996, Sect. 4.2)

$$\zeta(\theta) = R_0 [1 + h P_2(\cos \theta)], \quad (5.16)$$

where $\zeta$ is the altitude of the surface and $|h| \ll 1$. Such a form for the deformation ensures the conservation of volume at first order in $h$. Following Ragazzo and Ruiz (2015), we assume that the deformation is linear and we define the matrix $\mathcal{A} = \mathbb{I} + \mathcal{E}$ such that the vector position $\mathbf{r}$ of a point of $\mathbb{B}$ is displaced at the position $\mathcal{A}\mathbf{r}$ by the deformation. If $(\vec{e}_x, \vec{e}_y, \vec{e}_z)$ denotes the orthonormal cartesian basis, with $\vec{e}_z$ pointing towards the north pole of $\mathbb{B}$, then these vectors are clearly eigenvectors of $\mathcal{A}$. Hence, $\mathcal{A}$ is diagonal in the cartesian basis and looking at Eq. (5.16), we have

$$\mathcal{A}(\vec{e}_x, \vec{e}_y, \vec{e}_z) = \text{diag} \left( 1 - \frac{h}{2}, 1 - \frac{h}{2}, 1 + h \right). \quad (5.17)$$

The increment in potential energy due to the equatorial bulge reads (e.g. Wahr, 1996, Sect. 4.2)

$$V(r, \theta) = -\frac{3}{5} R_0^2 \omega_c^2 h \frac{R_0^3}{r^3} P_2(\cos \theta), \quad (5.18)$$
and, from Eqs. (5.12) & (5.13) and using $k_2(t) \propto \delta(t)$, the second Love number $k_2$ is

$$k_2(\eta) = -\frac{9h}{5\eta}. \qquad (5.19)$$

We need to determine the ratio $h/\iota$. A way to do so is to interpret $h$ as the generalized coordinate of the Lagrangian

$$\mathcal{L} = \mathcal{T}(h) - \mathcal{U}_g(h) - \mathcal{U}_{\text{el}}(h), \qquad (5.20)$$

and to obtain its value from the least action principle. In this Lagrangian, $\mathcal{T}$ is the kinetic energy of rotation of $\mathbb{B}$ while $\mathcal{U}_g$ and $\mathcal{U}_{\text{el}}$ are the gravitational binding energy and the elastic energy, respectively. For these quantities, we only consider the increment due to the rotation. A direct computation gives

$$\mathcal{T}(h) = \int_{0}^{2\pi} \int_{0}^{\pi} \int_{R_0}^{\zeta(\theta)} \frac{1}{2} \rho r^4 \omega^2 \sin^3 \theta dr d\theta d\varphi = -\frac{1}{5} m \omega_c^2 \iota h R_0^2, \qquad (5.21)$$

where $\rho$ is the density of $\mathbb{B}$. Denoting $\boldsymbol{\sigma}$ the stress tensor of $\mathbb{B}$ and $\boldsymbol{\varepsilon}$ its strain tensor, the elastic energy of deformation of $\mathbb{B}$ is

$$\mathcal{U}_{\text{el}} = \frac{1}{2} \int_{\mathbb{B}} \boldsymbol{\sigma} \cdot \boldsymbol{\varepsilon} \, d\boldsymbol{r} = \mu \int_{\mathbb{B}} \boldsymbol{\varepsilon} \cdot \boldsymbol{\varepsilon} \, d\boldsymbol{r}, \qquad (5.22)$$

where we used Hooke's law for a perfectly elastic body $\boldsymbol{\sigma} = 2\mu\boldsymbol{\varepsilon}$. The strain tensor is given by (Ragazzo and Ruiz, 2015)

$$\boldsymbol{\varepsilon} = \frac{1}{2} (\mathcal{A} + {}^t\mathcal{A}) - \mathbb{I} = \mathcal{E}. \qquad (5.23)$$

Equation (5.17) yields $\boldsymbol{\varepsilon} \cdot \boldsymbol{\varepsilon} = 3h^2/2$, and by conservation of the volume, we obtain

$$\mathcal{U}_{\text{el}}(h) = \frac{3}{2} \mu h^2 \int_{\mathbb{B}} dr = 2\pi \mu h^2 R_0^3. \qquad (5.24)$$

The gravitational binding energy is the potential energy of auto-gravitation

$$\mathcal{U}_g = -\frac{G\rho^2}{2} \int_{\mathbb{B}} \int_{\mathbb{B}} \frac{dr_1 dr_2}{|\mathbf{r}_1 - \mathbf{r}_2|} = -\frac{G\rho^2}{2} \int_{\mathbb{B}_0} \int_{\mathbb{B}_0} \frac{dr_1 dr_2}{|\mathcal{A}(\mathbf{r}_1 - \mathbf{r}_2)|}, \qquad (5.25)$$

where $\mathbb{B}_0$ is the ball of radius $R_0$ centered at $\mathcal{O}$. Ragazzo and Ruiz (2015) give

$$\mathcal{U}_g = -\frac{3m^2G}{10R_0} \int_0^{+\infty} \frac{d\lambda}{\sqrt{\det(\mathcal{A}^2 + \lambda\mathbb{I})}}. \qquad (5.26)$$

Given the simple form of $\mathcal{A}$ (Eq. (5.17)), this integral can be computed analytically, and at second order in $h$, we have

$$\int_0^{+\infty} \frac{d\lambda}{\sqrt{\det(\mathcal{A}^2 + \lambda\mathbb{I})}} = 2 - \frac{2}{5}h^2. \qquad (5.27)$$
The factor 2 gives the gravitational binding energy of a homogeneous ball and we discard it, since we are only interested in the increment due to the deformation. Thus we have

$$ U_g(h) = \frac{3}{25} m \omega_c^2 R_0^2 h^2. \qquad (5.28) $$

We normalize the Lagrangian $\mathcal{L}$ by the constant $-m\omega_c^2 R_0^2$ and obtain

$$ \mathcal{L} = \frac{1}{5} \iota h + \frac{3}{25} h^2 + \frac{3}{2} \frac{\mu}{\rho g R_0} h^2, \qquad (5.29) $$

where $g = Gm/R_0^2$ is the surface gravity of $\mathbb{B}$ when undeformed. The Euler-Lagrange equation (Eq. (2.12)) yields

$$ \frac{h}{\iota} = -\frac{5}{6} \left( 1 + \frac{25}{2} \frac{\mu}{\rho g R_0} \right)^{-1}, \qquad (5.30) $$

and then, using Eq. (5.19) (e.g. Roberts and Nimmo, 2008)

$$ k_2(\eta) = \frac{3}{2} \left( 1 + \frac{25}{2} \frac{\mu}{\rho g R_0} \right)^{-1}. \qquad (5.31) $$

In the original computation by Love, the coefficient 25/2 is 19/2 instead. The discrepancy is due to the assumption that the deformation is linear ($r \mapsto Ar$), which induces an overestimation of the elastic energy in this work (Ragazzo and Ruiz, 2015, Appendix 2).

For a perfectly fluid body ($\mu = 0$), we retrieve the known value $k_2(\eta) = 3/2$, or equivalently, $k_2(t) = 3\delta(t)/2$. The energy dissipated into heat due to tidal dissipation is proportional to the imaginary part of the second Love number, written in the frequency domain (Segatz et al., 1988, Eq. (13)). For a perfectly elastic body, $\mathfrak{Im}k_2 = 0$ implies the absence of energy dissipation. Such a behaviour is not satisfying because for realistic celestial bodies, the long term evolution of the orbits comes from the dissipation of the mechanical energy.

### 5.2.2 A simple viscoelastic body : the Maxwell model

Besides having a zero imaginary part, the second Love number in Eq. (5.31) was obtained using time-independent potentials (Eqs. (5.15) & (5.18)) and thus does not depend on the frequency $\eta$. Considering time dependent potentials, the form of Eq. (5.31) does not change (Appendix B of Efroimsky, 2012, Appendix A of Correia et al., 2014) and we have

$$ k_2(\eta) = \frac{3}{2} \left( 1 + \frac{25}{2} \frac{\mu(\eta)}{\rho g R_0} \right)^{-1}. \qquad (5.32) $$

The complex shear modulus $\mu(\eta)$ depends on the rheology of $\mathbb{B}$ and we consider, for the Maxwell model, a homogeneous body with a spring-like behaviour of shear modulus $\mu$ and a piston-like behaviour of damping coefficient $\tau$, put in serie (Christensen, 1982, Fig. 1.2).

<figure><img src="image_13.png" /></figure>
In the time domain, the differential equation verified by this system is $\sigma = 2\mu\epsilon + 2\tau\dot{\epsilon}$, while it reads $\sigma = 2\mu(1 + i\eta\tau)\epsilon$ in the frequency domain, where we defined the Maxwell relaxation time as $\tau = \tau/\mu$. Hence we have

$$ \mu(\eta) = \mu(1 + i\eta\tau). \qquad (5.33) $$

Other rheologies can be imagined by forming different mechanical systems. In the Kelvin model for instance, the spring and piston are in parallel and $\mu(\eta) = \mu i\eta\tau / (1 + i\eta\tau)$. The Maxwell model corresponds to a low-pass filter, since $\lim_{\eta \to +\infty} k_2 = 0$. In the time domain, its second Love number reads

$$ k_2(t) \propto e^{-t/\tau'} \quad \text{with} \quad \tau' = \tau \left( 1 + \frac{2 \rho g R_0}{25 \mu} \right)^{-1}, \qquad (5.34) $$

which corresponds to an exponentially decreasing memory of the past. Furthermore, $\mathfrak{Im}k_2(\eta) \neq 0$ ensures the dissipation of energy.

### 5.2.3 The constant-$\Delta t$ model

While the behaviour of the Maxwell model is qualitatively satisfying, the form of Eq. (5.34) makes it impossible to obtain an analytical expression of the perturbed potential $V$, because of the convolution product (5.12). In the limit when $\eta$ goes to zero, a first-order Taylor expansion in the vicinity of $\eta = 0$ in $k_2(\eta)$ (Eq. (5.32)) shows that there exists a scalar $\kappa_2$ and a time $\Delta t$, depending on the rheology of $\mathbb{B}$, such that the second Love number reads<fnref n="2" /> $k_2(\eta) = \kappa_2(1 - i\eta\Delta t) + \mathcal{O}(\eta^2\Delta t^2)$. The constant-$\Delta t$ model is thus defined as

$$ k_2(\eta) = \kappa_2 e^{-i\eta\Delta t} \approx \kappa_2 (1 - i\eta\Delta t), \qquad (5.35) $$

or equivalently

$$ k_2(t) = \kappa_2 \delta(t - \Delta t). \qquad (5.36) $$

The Dirac distribution in Eq. (5.36) solves the problem that the Maxwell model yields when trying to obtain an analytical expression for the perturbed potential. The constant-$\Delta t$ model is the limiting behaviour of any tidal model when $\eta \to 0$. In this model, the time $\Delta t$ corresponds to the memory of $\mathbb{B}$, that is, the time between the stress generated by the perturbing body and the response in the perturbed potential $V$. In the literature, the dimensionless parameter $\kappa_2$ is often referred to as the second Love number.

Due to the Taylor expansion, the constant-$\Delta t$ model is only valid when $\eta\Delta t \ll 1$. If $\eta\Delta t \gg 1$, $W$ evolves on timescales much shorter than $\Delta t$, and looking at Eq. (5.36), this means that $\mathbb{B}$ undergoes deformations on timescales much shorter than $\Delta t$. Since $\Delta t$ represents the time needed by $\mathbb{B}$ to reach hydrostatic equilibrium after a stress, the shape of $\mathbb{B}$ cannot evolve on timescales shorter than $\Delta t$. In that case, if $W$ were to suddenly start oscillating at high frequency, then $V$ would not react for a time $\Delta t$ and then would also suddenly start oscillating at high frequency, with an amplitude $\kappa_2$ that of $W$. This is obviously highly unrealistic. In practice, the constant-$\Delta t$ model is a good approximation as $\Delta t$ is at most of the order of a few minutes, while $2\pi/\eta$ generally exceeds one day.

<footnote n="2">²Any realistic tidal model is such that $\eta = 0 \Rightarrow \mathfrak{Im}k_2 = 0$, to ensure the absence of energy dissipation in the static case.</footnote>
### 5.2.4 Quality factor

The quality factor $Q$ is related to the amount of energy dissipated in a period $T = 2\pi/\eta$, and is defined as (MacDonald, 1964, Eq. (124))

$$Q^{-1} = \frac{1}{2\pi \mathfrak{E}^*} \int_{0}^{T} \frac{d\mathfrak{E}}{dt} dt, \qquad (5.37)$$

where $\mathfrak{E}(t)$ is the energy stored inside of $\mathbb{B}$ and $\mathfrak{E}^*$ is its maximum over the period $T$, called peak energy and defined as

$$\mathfrak{E}^* = -\mathfrak{E}(T/4) = -\int_{0}^{T/4} \frac{d\mathfrak{E}}{dt} dt. \qquad (5.38)$$

The stored energy $\mathfrak{E}(t)$ is given by the amount of work done by the perturbing potential $W$, that is

$$\frac{d\mathfrak{E}}{dt} = \int_{\mathbb{B}} \rho v \cdot \nabla W dr, \qquad (5.39)$$

where $v$ is the speed of the elemental volume surrounding $r$, in the frame attached to $\mathbb{B}$. The incompressibility of $\mathbb{B}$ reads $\nabla \cdot v = 0$ and yields $v \cdot \nabla W = \nabla \cdot (Wv)$. Using Green-Ostrogradsky theorem and neglecting the tangential contributions of $v$, Eq. (5.39) can be rewritten

$$\frac{d\mathfrak{E}}{dt} = \int_{\partial\mathbb{B}} \rho W(R_0, \sigma, t) \frac{\partial\zeta(\sigma, t)}{\partial t} d\sigma, \qquad (5.40)$$

where $\partial\mathbb{B}$ is the surface of $\mathbb{B}$, $\zeta$ its altitude (Eq. (5.16)), $\sigma = (\theta, \varphi)$ and $d\sigma = \sin\theta d\varphi d\theta$. If $\mathbb{B}$ is excited with a frequency $\eta$, that is

$$W(r, \sigma, t) = W_0(r, \sigma) \cos \eta t, \qquad (5.41)$$

then, there exists a phase lag $\delta$, accounting for the delay in the response due to the non-elasticity of $\mathbb{B}$, such that the altitude of the surface takes the form

$$\zeta(\sigma, t) = \zeta_0(\sigma) \cos(\eta t - \delta). \qquad (5.42)$$

Equation (5.40) yields $d\mathfrak{E}/dt = -\eta\rho\mathfrak{U} \cos\eta t \sin(\eta t - \delta)$, where $\mathfrak{U} = \int_{\partial\mathbb{B}} W_0(R_0, \sigma)\zeta_0(\sigma)d\sigma$. The energy dissipated over one period and the peak energy then read

$$\int_{0}^{T} \frac{d\mathfrak{E}}{dt} dt = \pi \rho \mathfrak{U} \sin \delta \quad \text{and} \quad \mathfrak{E}^* = \frac{\rho \mathfrak{U}}{2} \cos \delta. \qquad (5.43)$$

Posing $\delta = \eta\Delta t$ and using Eq. (5.37), we finally have for $Q$ (MacDonald, 1964, Eq. (130))

$$Q^{-1} = \tan \delta = \tan \eta \Delta t \approx \eta \Delta t. \qquad (5.44)$$

In the literature, it is often assumed that $Q$ is constant, instead of $\Delta t$ (e.g. Lainey et al., 2012, Fig. 2). In the rest of this manuscript, we adopt the constant-$Q$ model<fnref n="3" /> and make use of the parameters $\kappa_2$ and $Q$. The timescales of orbital decay due to tides do not depend on both $\kappa_2$ and $Q$, but rather on the ratio $\kappa_2/Q$ (see Eq. (6.29)). We give this ratio in Table 5.1 for some bodies of the Solar System.

<footnote n="3">³With a unique frequency of excitation $\eta$, the constant-$\Delta t$ and constant-$Q$ models are identical.</footnote>