Issue 
A&A
Volume 664, August 2022



Article Number  A1  
Number of page(s)  19  
Section  Celestial mechanics and astrometry  
DOI  https://doi.org/10.1051/00046361/202243261  
Published online  04 August 2022 
Dynamics of coorbital exoplanets in a firstorder resonance chain with tidal dissipation
^{1}
IMCCE, UMR8028 CNRS, Observatoire de Paris, PSL Univ., Sorbonne Univ.,
77 av. DenfertRochereau,
75014
Paris, France
email: jeremy.couturier@obspm.fr
^{2}
CFisUC, Departamento de Fisica, Universidade de Coimbra,
3004516
Coimbra, Portugal
Received:
3
February
2022
Accepted:
8
April
2022
Coorbital planets (in a 1: 1 mean motion resonance) can be formed within a Laplace resonance chain. We develop a secular model tc study the dynamics of the resonance chain p: p : p + 1, where the coorbital pair is in a firstorder mean motion resonance with the outermost third planet. Our model takes into account tidal dissipation through the use of a Hamiltonian version of the constant timelag model, which extends the Hamiltonian formalism of the pointmass case. We show the existence of several families of equilibria, anc how these equilibria extend to the complete system. In one family, which we call the main branch, a secular resonance between the libration frequency of the coorbitals and the precession frequency of the pericentres has unexpected dynamical consequences when tidal dissipation is added. We report the existence of two distinct mechanisms that make coorbital planets much more stable within the p : p : p + 1 resonance chain rather than outside it. The first is due to negative real parts of the eigenvalues of the linearised system with tides, in the region of the secular resonance mentioned above. The second comes from nonlinear contributions of the vector fielt and is due to eccentricity damping. These two stabilising mechanisms increase the chances of a future detection of exoplanets in the coorbital configuration.
Key words: celestial mechanics
© J. Couturier et al. 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model. Subscribe to A&A to support open access publication.
1 Introduction
Coorbital systems have been extensively studied since the discovery of equilibria in the threebody problem by Euler (1764) and Lagrange (1772). In the hierarchical case, that is when (m_{1} + m_{2}) /m_{0} < 1/27, where m_{0} is the mass of the central body and m_{1} and m_{2} are the coorbital masses, Gascheau (1843) showed that the equilateral equilibria, where the bodies are at the vertices of an equilateral triangle, are linearly stable. This result, combined with the discovery of several examples of coorbital bodies, such as the Jovian trojans or the horseshoeshaped orbits of Janus and Epimetheus around Saturn, contributed to increasing the interest of scientists for systems of this kind.
In the planar and circular case, for masses in the range 3 × 10^{−4} < (m_{1} + m_{2}) /m_{0} < 1/27, the angle λ_{1} – λ_{2} between the coorbitals librates around its equilibrium of ±60° in trajectories commonly called tadpole orbits, where the coorbital angle is bounded by^{1} 23.9° < λ_{1} – λ_{2} < 180°. However, for masses (m_{1} + m_{2}) /m_{0} < 3 × 10^{−4} (e.g. Laughlin & Chambers 2002), a separatrix in the phase space delimits a region of stable trajectories of another kind, generally called horseshoe orbits, where the critical angle λ_{1} – λ_{2} librates around 180° with at least 312.2° of amplitude. For low eccentricities and small libration amplitudes, still in the planar case, the existence of two proper modes called Lagrange and antiLagrange have been shown, numerically by Giuppone et al. (2010) and analytically by Robutel & Pousse (2013). In the Lagrange configuration, the pericentres of the orbits do not precess and verify the relation ϖ_{1} – ϖ_{2} = 60°, whereas in the antiLagrange configuration, both orbits precess at the same frequency while maintaining the relation ϖ_{1} – ϖ_{2} = 240^{°}. For low eccentricities but very large libration amplitudes, in the region of horseshoeshaped orbits, Couturier et al. (2021) showed that the Lagrange and antiLagrange configurations correspond to ϖ_{1} – ϖ_{2}= 0° and ϖ_{1} – ϖ_{2} = 180°, respectively. High eccentricities give rise to topological changes in the phase space (see Leleu et al. 2018), and thus to many more exotic trajectories, while the dynamics of the inclined problem is even more complex by allowing, among other things, transitions between these orbits and the retrograde coorbitals (Namouni 1999).
The discovery of exoplanets raised the question of the existence of coorbital planets, which are absent from the Solar System. Accretion in situ at the equilateral equilibria of a primary or capture in the 1:1 resonance of planets formed in other parts of the system are two possible scenarii of formation of such systems (e.g. Laughlin & Chambers 2002; Cresswell & Nelson 2009). The stability of coorbital planets formed in a disc has been studied by Leleu et al. (2019), who show that under dissipative interactions with the gas disc, the equilateral equilibria can be either attractive or repulsive, depending on the coorbital mass ratio and the parameters of the disc. Moreover, Leleu et al. (2019) show that, at least around lowmass stars, coorbital exoplanets generally end up in a tadpole configuration and often within a Laplace resonance chain.
For coorbital exoplanets orbiting close to the host star, tidal dissipation induced by the differential gravitational interaction leads to a longterm evolution of the orbits. Couturier et al. (2021) shows that, for a pair of coorbital exoplanets orbiting a star, the equilateral Lagrangian equilibria are always repulsive under tidal interactions, and that regardless of the parameters of the system, the destruction of the coorbital motion is unavoidable. However, the discovery of coorbital exoplanets is still possible because the destruction time is strongly dependent on the parameters of the system and can easily be greater than the lifetime of the host star. Couturier et al. (2021) neglected any interaction with possible other planets in the system. In this paper we extend this work to the case where the pair of coorbital exoplanets interacts with an outermost third planet, in a firstorder mean motion resonance with the coorbitals. More precisely, we study the Laplace resonance chain p : p : p + 1, where p is a small integer.
In Sect. 2, we study the pointmass p : p : p + 1 resonance chain in the absence of tides. We show how rich and complex the dynamics of this chain is, and we conclude the section with the presentation of the stability map of the chain. In Sect. 3, we include tidal dissipation in the model by an extension of the Hamiltonian formalism. We study the linearised system in the vicinity of the equilibria, and by computing the real parts of the eigenvalues, we show the existence of a linearly stable zone around the 1 : 1 secular resonance between the libration frequency of the coorbital and the precession frequency of the pericentres. In Sect. 4, we compare the analytical results with numerical simulations. They confirm the results of Sects. 2 and 3 and highlight the existence of a stabilisation mechanism of the coorbitals due to eccentricity damping. We discuss our results in Sect. 5. In Table A.1, we list the notations used throughout this paper. Appendix F completes Sect. 4 with more numerical simulations and a complete discussion on the influence of the mass of the third planet on the coorbital dynamics.
2 . The p : p : p +1 Resonance Chain
2.1 The Hamiltonian of the Problem
In this section, we study an occurrence of the pointmass planar fourbody problem. We construct the Hamiltonian associated with the resonance chain p : p : p + 1, where a central body, a star of mass m_{0}, is orbited by two coorbital planets of mass m_{1} and m_{2} and a third planet of mass m_{3} further away from the star and in a firstorder mean motion resonance with the pair of coorbital planets. Although we write all equations for a general value of the integer p, the figures are restricted to the case p = 1 where the nominal period of the third planet is twice that of the coorbitals. For all planets we define the quantities β_{j} = m_{0}m_{j}/(m_{0} + m_{j}) and , where is the gravitational constant.
2.1.1 The Averaged Hamiltonian
In order to define a canonical coordinate system related to the semimajor axis a_{j}, the eccentricity e_{j}, the mean longitude λ_{j}, and the longitude of the pericentre ϖ_{j} of planet j, we first consider Poincaré heliocentric coordinates , where (1)
In these coordinates, the Hamiltonian derives from the symplectic form (2)
Following Laskar & Robutel (1995), the planetary Hamiltonian is written (3)
where the Keplerian part, due to starplanet interactions, reads (4)
and the perturbation H_{P}, whose size relative to H_{K} is of order (m_{1} + m_{2} + m_{3}) /m_{0}, is expanded in power series of the eccentricities. We assume that the system is close to the resonance p : p : p + 1. This means that the nominal mean motions verify (5)
while the nominal semimajor axes a_{j,0} are related to n_{j,0} by the Kepler law . The a_{j} are always close to their nominal value^{2} a_{j,0}, and the stay close to the quantities defined as (6)
To study the dynamics in the vicinity of the resonance, we expand the Hamiltonian in the neighbourhood of . An expansion at order 2 in the Keplerian part and at order 0 in the perturbative part generates remainders of the same size, and we limit ourselves to (7)
A suitable linear change of variables to deal with the p : p : p + 1 resonance chain is (e.g. Delisle 2017) (8)
which is canonical if we transform the actions according to (9)
Since the total angular momentum G′ is a first integral, the Hamiltonian does not depend on the angle ξ_{3}. Moreover, in the p : p : p + 1 resonance, the angle ξ_{2} is fast circulating and we average over it. The averaged Hamiltonian reads (10)
This change of variable, along with the averaging process, allows the two degrees of freedom associated with (ξ_{2}, ξ_{3}, Γ′, G′) to be lost, and we are left with four degrees of freedom. After the averaging process, the scaling factor Γ′ and the angular momentum G′ are both parameters and a rescaling by Γ′ reduces the dependency to only one parameter. As we study the effect of tidal dissipation on the dynamics in Sect. 3, it is actually more convenient to normalise by the constant rather than by Γ′, which is not constant when dissipation is present. In other words, we perform the canonical transformation (11)
while the angles are unchanged.
2.1.2 Expansion of the Keplerian Part
As mentioned above, the Keplerian part of the Hamiltonian is expanded at second order in the vicinity of the . If we denote , the expansion reads (12)
Substituting for , we obtain (13)
The third term is constant and can be removed without changing the dynamics. Performing the normalisation (Eq. (11)) and the change of variable (Eq. (9)), we have (14)
where we denoted Υ = G + D_{1} + D_{2} + D_{3} = ∑_{j}Λ_{j} and , that is (15)
Without dissipation, Γ is constant and almost equal to 1 and we simply evaluate ℋ_{k} at Γ = 1, hence achieving the reduction to only one parameter^{3}. In this case, the last term is constant and can also be removed.
Instead of the variables L, G, and Γ, we can use the variables ΔL, ΔG, and ΔΓ defined by their difference to the Keplerian resonance (Eq. (6)). In this case, the approximation Γ = 1 becomes ΔΓ = 0 and Eq. (12) yields, once normalised, (16)
where ΔΥ = ΔG + D_{1} + D_{2} + D_{3} = ∑_{j} ΔΛ_{j}. We find the Hamiltonian (16) to be well adapted to the analytical work derived in Sect. 2.2.2, while we use the Hamiltonian (14) in the remaining sections. Moreover, we do not perform the evaluation Γ = 1 in Sect. 3, where tidal dissipation is present and Γ is a variable quantity. Both Hamiltonians yield the same dynamics and it is only a matter of preference.
2.1.3 Expansion of the Perturbative Part
The perturbative part ℋ_{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 (17)
We denote For a pair (p_{1}, p_{2}) ϵ {(1,2), (1,3), (2,3)} of planets, the perturbation to the Hamiltonian due to their mutual interaction reads (Laskar & Robutel 1995) (18)
For a nonzero , the conservation of the angular momentum imposes on the tuples q = (q_{1}, q_{2}, q_{3}, q_{4}) ϵ ℕ^{4} and k = (k_{1}, k_{2}) ϵ ℤ^{2} to verify the d'Alembert rule: (19)
This rule, combined with the averaging process, implies that ℋ_{1,2} has no odd term in eccentricity, while ℋ_{1,3} and ℋ_{2,3} have no term of order 0. Since we limit ourselves to the second order in eccentricity, we write (20)
where the superscript refers to the order in eccentricity, while the subscript refers to the considered pair of planets. The Hamiltonian ℋ^{(0)} has no subscript since only the pair of coorbitals yields terms of order 0 and no confusion is possible.
Following Laskar & Robutel (1995), ℋ_{1,3}and ℋ_{2,3}can be written (21)
where = 1 if j = 1 and zero otherwise. The quantities depend only on p and can be obtained using the Laplace coefficients. For p = 1, their analytical expressions, as well as a numerical evaluation, is given in Appendix B.
The perturbation ℋ_{1,2}cannot be obtained using the same procedure since the Laplace coefficients diverge for equal nominal semimajor axes and a_{1,0} = a_{2,0} = ā. We instead follow the method described in Robutel & Pousse (2013). We denote and find (23)
The final simplified Hamiltonian is then (25)
We denote F_{0} : ℝ^{8} ↦ ℝ^{8} the differential system derived from Eq. (25) by the HamiltonJacobi equations.
2.2 Equilibria and Linearisation in Their Vicinity
In this section, we study the equilibria of the resonance p : p : p + 1 and the dynamics in their vicinity.
2.2.1 Fixed Points and Libration Centres
One of the consequences of averaging over the mean motion is that the averaged Hamiltonian (25) has equilibria, that is, points in the phase space where its gradient vanishes. The complete Hamiltonian (3) has no equilibria, however, and since the Hamiltonian (25) is supposed to model it, we are interested in the dynamics of the complete Hamiltonian at the equilibria of the model.
At a fixed point (or equilibrium) of the model, L and D_{j} are constant and so are e_{j} and a_{j}. Similarly, the angles σ_{j} and ξ are constant; that is, there are constants c_{j} such that (26)
However, the secular angle –pλ_{2} + (p + 1) λ_{3} and the pericentres a j are not constant at the equilibria, but they all precess with the same frequency, which we denote ν_{3}.
The average performed in Eq. (10) is actually analogous to a firstorder Lie serie expansion, and the averaging process can be seen as a periodic change of variable. Denoting x and x′ the variables of the Hamiltonian respectively before and after the average, we have , where L_{W} = {W, ·} denotes the total time derivative along the trajectories of the scalar field W, which is constrained by the cohomological equation (Deprit 1969) (27)
This equation shows that, at the equilibria, W is periodic of time, and so is the change of variable. This means that fixed points in the averaged model correspond to periodic trajectories in the complete system. For a quantity not invariant by rotation around the axis of the total angular momentum, however, like the secular angle –pλ_{2} + (p + 1) λ_{3} or the pericentres ϖ_{j}, a fixed point in the model corresponds in the complete Hamiltonian to a quasiperiodic motion with the two frequencies ν_{2} and ν_{3}, with (28)
where these qua tities are evaluated at the equilibrium. More precisely, it is a periodic motio with frequecy ν_{2} i a rotatig frame followig all the perice tres at frequecy ν_{3}. This result holds true for any resoace chain (e.g. Eq. (A.1) in Delisle 2017). In the rest of this work, what is known as a fixed point or equilibrium for the model is referred to as a libratio cetre in the complete system.
2.2.2 Analytical Results at First Order in Eccentricity
Even if truncated at order 1 in eccentricity, the fixed points of the Hamiltonian (25) cannot be given analytically. Similar difficulties were met by Delisle (2017) for resonance chains with firstorder resonances between nonconsecutive planets. However, we show here that a further simplification allows us to obtain analytical expressions of the equilibria and of the eigenvalues of the linearised system.
To further simplify the Hamiltonian, we force a decoupling between the degree of freedom (ξ, ΔL) associated with the libration of the coorbitals and the three other degrees of freedom (σ_{j}, D_{j}) To this end, we first evaluate at ξ = π/3. Indeed, the Hamiltonian ℋ_{K} + ℋ^{(0)} + ℋ^{(1)} is only a perturbation of ℋ_{K} + ℋ^{(0)} that only has equilibria at ξ ϵ {±π/3, π}, where ξ = π is the hyperbolic (unstable) aligned configuration and ξ = ±π/3 are the equilateral elliptic (stable) equilibria (Robutel & Pousse 2013). Both elliptic equilibria are symmetric with the same dynamics, hence we only consider ξ = π/3. Then, we replace the variable ΔL by the constant ΔL^{★} in the antidiagonal term ^{4} of ^{ℋ}_{K} in (16), where (29)
is the value of ΔL for which ∂ℋ_{K}/∂ΔL vanishes. The D_{j,0} are given by Eq. (32). While the evaluation at ξ = π/3 allows analytical expressions for the position of the fixed points, the evaluation at ΔL = ΔL^{★} also uncouples (ξ, ΔL) from (σ_{j}, D_{j}) and enables analytical expressions of the eigenvalues of the linearised system in the vicinity of the fixed points. The differential system derived from ℋ_{K} + ℋ^{(0)} + ℋ^{(1)} once these simplifications have been performed, is given in Appendix C. It vanishes when the angles are equal to (30)
and when the actions are equal to^{5} (32)
where the precession frequency of the pericentres is (33)
and we define the constant H by (34)
The unknowns of Eq. (32) are the Dj,0, and since the ratios D_{2,0}/D_{1,0} and D_{3,0}/D_{1,0} are known, we are reduced to the unique unknown D_{1,0}. The precession frequency of the pericentre, ν_{3}, depends on D_{1,0} (see Eqs. (29) and (33)). Denoting C = 1 + D_{2,0}/D_{1,0} + D_{3,0}/D_{1,0} and performing the translation Z = D_{1,0} + 2ΔG/3C, Eq. (32) is rewritten as a thirddegree polynomial in Z:
Z^{3} – PZ – Q = 0, where (35)
The coefficients P and Q of this polynomial depend on the parameter ΔG. There is a bifurcation between 1 and 3 real solutions when 27Q^{2} – 4P^{3} = 0, that is at (36)
In the rest of this work, we use the parameter (37)
In this section, Γ = Γ′/Γ^{★} ≈ 1 is simply evaluated at 1 and ignored, but not in Sects. 3 and 4, where tidal dissipation induces a drift in Γ, hence in δ. The normalisation by ΔGbif ensures that the bifurcation is at δ = 1, regardless of the planetary masses.
The forced decoupling that we performed to obtain these expressions leads to results that are very similar to the second fundamental model of resonance proposed by Henrard & Lemaitre (1983). The fixed points are given by the roots of the thirddegree polynomial in Z (Eq. (35)), which has one or three real solutions depending on δ, hence a bifurcation. The solutions of Eq. (35) are plotted in Fig. 1. For δ < 1, only one elliptic equilibrium exists, called the main branch, while for δ ≥ 1 two other fixed points appear, one of them being hyperbolic, hence the presence of separatrices in the phase space and the formal existence of a resonance. These results come from a strong hypothesis, and we see in Sect. 2.2.3 that the topology of the Hamiltonian (25) is different (see Table 1). However, we show in Fig. 2 that these analytical expressions are accurate for low eccentricities.
In the vicinity of the main branch, we linearise the differential system. We use the cartesian coordinates (39)
and denoting X = ^{t}(u_{1},u_{2},u_{3},υ_{1},υ_{2},υ_{3}, ΔL, ξ), the linearised systemis (40)
and X_{0} is the equilibrium value of X. The matrix Q_{6} is given in Appendix D. Its characteristic polynomial reads (41)
where I is defined in Eq. (D.1). It is interesting to note that the precession frequency of the pericentres, ±iv_{3}, is an eigenvalue of the differential system (Eq. (40)). This factorisation was already noted by Pucacco (2021), who studied the resonance chain 1 : 2 : 4 of the Galilean satellites, although it was not attributed to the precession of the pericentres. The eigenvalues of Q_{2} are ±iv, where (42)
is the libration frequency of the angle ξ in the neighbourhood of the equilateral Lagrangian configuration (Robutel & Pousse 2013; Couturier et al. 2021). Figures 1 show that ν3 < 0 for the main branch, which ensures that the roots of Eq. (41) are purely imaginary. Evaluating the eigenvalues ±iv and ±iv_{3} in the vicinity of the main branch shows that, at δ ≈ −5.6 for the planetary masses in Fig. 1, all these eigenvalues have roughly the same value, yielding a 1 : 1 secular resonance between the libration frequency of the coorbitals and the precession frequency of the pericentres. Other secular resonances between ν and ν3 are shown in Fig. 1, and are also very visible on the stability map from Fig. 6. We show in Sect. 3.2 that the secular resonance 1 : 1 has important consequences for the tidal stability of the coorbital pair.
Fig.1 Values of the eccentricities and fundamental frequencies given by the analytical results. Left: value of e_{1} at the fixed points of the resonance chain 1:1:2 as predicted by Eq. (35). Right: values of ν/η and ν_{3}/η along the main branch as a function of δ, for the same resonance chain, predicted by Eqs. (42) and (33). In both panels the planetary masses are (m_{1} + m_{2}) /2 = m_{3} = 10^{−4} m_{0} and m_{2}/m_{1} = 10. According to Eq. (32), e_{1} = e_{2}, and for this choice of masses and resonance chain e_{3}/e_{1} = 0.5471. In Table 1, which gathers the fixed points at second order in eccentricity, the three branches visible in the left plot are branches 1 (main branch), 2, and 6. The secular resonances between ν and ν_{3} are shown on the right, and are also easy to spot in Fig. 6. 
The 15 equilibria of the simplified Hamiltonian (25).
2.2.3 Topology of the Phase Space
Limiting the work at first order in eccentricity and forcing a decoupling between (L, ξ) and (D_{j}, σ_{j}) gives analytical expressions of the linearised system^{6}, but at the cost of strong approximations. We develop here a Newton–Raphson–based algorithm to numerically find the equilibria of the model (25) without these approximations.
The vector field F_{0} derived from the Hamiltonian (25) depends on the choice of the parameter δ (through ΔG), and once a fixed point is found for a particular value of δ, we repeat the NewtonRaphson algorithm for slowly varying values of δ in order to travel along the whole branch. We look for equilibria exploring the parallelepiped in the phase space defined byu_{j} < 0.08 and . We choose as initial condition of the NewtonRaphson method for L since all equilibria are expected to be close to this value, and so no discretisation is necessary along this axis. In the same way, we only choose ξ_{0} ϵ {±π/3, π}.
We display in Table 1 all the equilibria that we have found for δ = 7, their hyperbolic or elliptic nature, and the value of δ that gives birth to the branch. Due to the difficulty of exploring a thinly discretised grid in six dimensions, we may not have found all the possible equilibria. We discretised the axes u_{j} and υ_{j} with only eight points, testing 3 × 8^{6} initial conditions, all of which converged towards 15 equilibria. Since the Hamiltonian (25) is invariant by the transformation (ξ, σ_{j}) ↦ (–ξ, –σ_{j}), fixed points with values of the angles different from 0 or π have a symmetric, hence the ± and ∓ signs in Table 1 (the upper sign corresponds to a fixed point and the lower sign to its symmetric). This symmetry corresponds to the invariance of the system by a rotation of angle π around an axis normal to the total angular momentum.
As can be seen in Table 1, only the branches 1, 2, and 3 of fixed points can be elliptic, and thus we focus only on them in the rest of this work. In Fig. 2 we plot these branches for values of δ ranging from −2 to 9. For branches 1 and 2, which are predicted by the first order in eccentricity, for comparison we also plot them as given by Eqs. (30) and (35). This section shows how the analytical model is unable to locate the equilibria of the simplified Hamiltonian (25) for e_{j} ≳ 0.005 (see Fig. 2) and does not even give its topology for e_{j} ≳ 0.05 (branches 3, 5, and 7 do not exist at first order in eccentricity; see Table 1). This discrepancy between first and second order in eccentricity was already mentioned by Beaugé et al. (2006) in the case of the twoplanet 1 : 2 mean motion resonance.
Fig. 2 Position of the elliptic branches 1, 2, and 3 (see Table 1) of equilibria of Hamiltonian (25) in the resonance 1 : 1 : 2, for −2 ≤ δ ≤ 9. For each branch three curves appear, corresponding to for j = 1,2,3. The planetary masses are the same as in Table 1. A zoomedin image close to the origin is shown in the inset. In this area, the analytical position of these equilibria, given by Eqs. (30) and (35), is plotted by a thin grey line. The agreement is good at low values of eccentricity, but quickly worsens further from the origin. In particular, the thin grey lines are straight since σ_{j}· does not depend on δ in Eq. (30). Branch 1 exists for all values of δ and has all colours from yellow to dark purple, while branch 3 only exists at δ > 5.997 and thus only has purple. 
2.2.4 Comparison with the Complete Hamiltonian
In this section we compare the positions of the equilibria of the secular (simplified) Hamiltonian (25) to those of the corresponding periodic orbits of the complete (full) Hamiltonian (3), which we call libration centres in Sect. 2.2.1. To this end, we develop an iterative algorithm, similar to that in Couetdic et al. (2010), to find a libration centre of the complete Hamiltonian using an equilibrium of the simplified Hamiltonian as initial condition.
We assume, close enough to a libration centre of the complete Hamiltonian (3), that the trajectories are quasiperiodic, and we write, for any complex quantity z depending on these trajectories, (43)
where the coordinates of ω = ^{t}(ν, ν_{2}, ν_{3}, g_{1}, g_{2}, g_{3}) are the fundamental frequencies of the complete Hamiltonian (3). The frequencies ν, ν_{2}, and ν_{3} have approximate values given by Eqs. (42), (28), and (33), respectively.
As explained in Sect. 2.2.1, the libration centres correspond to points in the phase space where the motion is periodic^{7} with frequency ν_{2}. The description of the algorithm is as follows:
For a choice of the parameter δ, find the position of an elliptic equilibrium of the simplified Hamiltonian (25) with a Newton–Raphson method. Use it as the initial condition to integrate numerically the trajectories of the complete Hamiltonian (3);
For a complex quantity z depending on the trajectories of Eq. (3), obtain the decomposition (Eq. (43)) using a frequency analysis method (e.g. Laskar 1993);
Identify the terms depending on frequencies other than ν_{2} and set to 0 the corresponding coefficient z_{k}.
Proceed similarly for different quantities z and evaluate them at time t = 0 in order to obtain a new initial condition. Restart from the first step using the new initial condition instead of the equilibrium of Eq. (25).
The process is iterated until a convergence occurs. In Fig. 3, we display the trajectories of the quantities in the plane (e_{j} cos σ_{j}, e_{j} sin σ_{j}) as the algorithm iterates. Isolating terms featuring frequencies other than ν_{2} is difficult, if not impossible, since the frequency analysis gives the scalars k•ω but not the vector k, which cannot be deduced as the vector ω is unknown for the complete Hamiltonian. We can get around this difficulty because ν_{2} is much larger than the other frequencies, hence it is easy to isolate terms that depend on ν_{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 ν_{2}.
This algorithm is only able to find libration centres associated with elliptic fixed points. We use it to find the branches of libration centres associated with branches 1, 2, and 3 of Table 1. We confirm the existence of a small hyperbolic zone for branch 2 in the complete Hamiltonian when the algorithm stops converging as it travels along the branch (by slowly incrementing the value of δ). We plot in Fig. 4 the main branch (branch 1) of the libration centres of the complete Hamiltonian (3) and we compare it with the main branch of equilibria of the simplified Hamiltonian (25). In Fig. 5, we plot the precession frequency of the pericentres, ν_{3}, for the main branch of the libration centres, and we compare it with the analytical expression (33). This figure shows that, for high values of δ, n_{1}/n_{3} tends towards (p + 1) /p (which is equal to 2 in this case). For small values of δ though, this ratio diverges from its nominal value and libration centres at small values of δ are far from the Keplerian resonance. For branch 1, we now refer to very negative values of δ as far from the resonance and to very positive values of δ as deep in the resonance. The value of ν_{3} in the complete Hamiltonian is obtained from the frequency analysis of , once the libration centre is known.
Fig. 3 Trajectories of in the complete Hamiltonian (3) for the six iterations needed for the algorithm to converge to the libration centre. Iteration 0 is the equilibrium (main branch) of the simplified Hamiltonian (25) at δ = 5. It is rather far from the libration centre of the complete Hamiltonian (see also Fig. 4). After six iterations, the algorithm converges to the libration centre and the motion is periodic (hence the closed curves) with frequency ν_{2} (see Sect. 2.2.1). The planetary masses are the same as in Table 1 and the resonance chain is 1:1:2. 
2.3 Stability Map of the p : p : p + 1 Resonance Chain
Before taking into account tidal dissipation in the model, we study the stability of the pointmass p : p : p + 1 resonance chain by constructing a dynamical map using the frequency analysis method to determine the chaoticity of a given orbit (Laskar 1990). More precisely, we study the stability of the chain along its main branch (see Fig. 1) between δ = 7, deep in the resonance, and δ = −7, outside the resonance.
For each value of δ, we compute the position of the exact libration centre by means of the algorithm described in Sect. 2.2.4, and we choose an initial value for the angle ξ between its equilibrium value (near 60°) and the value at the boundary between tadpole and horseshoeshaped orbits (close to 24°, see Robutel & Pousse 2013). For values of ξ_{0} close to 60°, the considered orbit is close to the main branch and it moves away for decreasing values of ξ_{0}. For all other variables, we choose as the initial condition the value at the libration centre. Every trajectory (i.e. every choice of δ and ξ_{0}) is integrated over 80 000 periods of the coorbital planets, and for each half of the simulation the exact value of the libration frequency ν^{(2)} is extracted from the frequency analysis of e^{iξ}. We obtain two values of ν, namely ν^{(1)} for the first half of 40 000 periods and ν^{(2)} for the second half. The diffusion index, defined as (Robutel & Gabern 2006) (44)
measures the degree of quasiperiodicity of the orbit. Orbits with ζ_{ν} < −6 are considered close to quasiperiodic (stable), while orbits with ζ_{ν} > −2 are very chaotic (unstable). We plot the stability map for the resonance chain 1 : 1 : 2 in Fig. 6. Secular resonances between ν and ν_{3}, already predicted by the analytical results in Fig. 1 are visible and induce chaotic motion. Overall, this stability map shows that the chain p : p : p + 1 is mainly stable.
Fig. 4 Position of branch 1 of elliptic libration centres of the complete Hamiltonian (3) in the resonance 1:1:2 for −7 ≤ δ ≤ 7. The planetary masses are the same as in Table 1. As a comparison, branch 1 of equilibria of the simplified Hamiltonian (25) is plotted in grey for the same range in δ. For small enough values of δ the eccentricity is not too high and the agreement is good. The simplified and the complete Hamiltonian diverge when δ → +∞. 
Fig. 5 ν_{3}/η as a function of δ, in the complete Hamiltonian, for the resonance chain 1 : 1 : 2, along branch 1 of the libration centres in Fig. 4. The analytical expression (33) is plotted in grey for comparison. The colour gives the value of n_{1} /n_{3}. As expected from Eq. (33), ν_{3} is proportional to δ far from the resonance where the eccentricities are small. 
3 Tides in the p : p : p + 1 Resonance Chain
In Sect. 2, we assumed that the bodies are point mass objects. In this section the approximation is removed and tidal dissipation due to differential and inelastic deformations of the bodies is taken into account.
3.1 Extended Hamiltonian and Equations of Motion
Tidal contributions to the orbital evolution of the system follow a very general formulation initiated by Darwin (1880). Differential interactions between the bodies raise tidal bulges and the subsequent redistribution of mass is responsible for a perturbation in the gravitational potential generated by body j at any point r in the space. This perturbation is given by (e.g. Kaula 1964) (45)
where the indice i (resp. j) refers to the body responsible for the tidal bulge (resp. where the bulge is raised), r_{i} is the position of body i with respect to the barycentre of body j, R_{j} is the radius of the body is its second Love number, P_{2} is the second Legendre polynomial, and S is the angle between and r. For a body of mass m_{k}, located at r_{k} and interacting with this bulge, the increment in potential energy is (46)
For N = 4 tidally interacting bodies, N (N − 1)^{2} = 36 such potentials are generated. In the case of planets orbiting a solartype star, however, only tides raised by the star on the planets and felt by the star should be considered since they are dominant with respect to any other contribution (see Couturier et al. 2021), which means that we only consider the three contributions U_{i,j,k} with i = k = 0 and j ∊ {1, 2, 3}.
The dissipation of mechanical energy inside the planets introduces a time delay ∆t between the tidal stress and the corresponding deformation. As a consequence, the tidal bulge is not aligned with the star and the subsequent torque affects the spins and orbits of the planets. For any quantity z_{j} related to planet j, we write (47)
For a frequency of excitation , the quality factor (Munk & MacDonald 1960), which measures the amount of energy dissipated in a period , is related to the time delay as (48)
The dependency of on is unknown, and a simple commonly used rheology consists in considering that the time delay is independent of the frequency (Mignard 1979). We adopt this tidal model in this work, reducing the rheology to the constant parameters and .
Although tides do not preserve the total energy, the Hamiltonian formalism is extended by considering the starred variables as parameters when deriving the equations of motion. Their contribution to the Hamiltonian reads (49)
where in the heliocentric reference frame (50)
and θ_{j} is the rotation angle of body is its rotation rate, is the conjugated momentum of θ_{j}, and α_{j} is a dimensionless structure constant depending on the state equation of body j such that is its principal moment of inertia. The transformations (9) and (11) are performed on the tidal Hamiltonian, with the normalisations . Denoting (52)
where ϑ stands for any angle, we obtain for the tidal Hamiltonian (53)
with (Couturier et al. 2021) (54)
We note that no expansion at order 0 in the vicinity of is performed since it loses relevant tidal dynamics (see Couturier et al. 2021). We instead keep exact expressions in Λ_{j}. The differential system is derived from the tidal Hamiltonian using the HamiltonJacobi equations and considering the starred variables as parameters. The starred variables are then expressed by a firstorder Taylor expansion in Eq. (47) (see Couturier et al. 2021 for more details). The perturbation to the vector field, due to tides and at second order in eccentricity, reads (56)
where we posed Since the differential system does not depend on ξ_{2} and ξ_{3} and as their dynamics are of no interest to us, the lines ξ_{2} and ξ_{3} are absent from the differential system (56).
Fig. 6 Diffusion index ζ_{ν} as a function of δ and ξ_{0}. The planetary masses are as in Table 1 and the chain is 1:1:2. The top of the figure, at ξ_{0} ≈ 60°, is the main branch. The bluetogreen regions are almost quasiperiodic (stable), while the red regions are chaotic (unstable). The secular resonances between the libration frequency ν and the precession frequency of the pericentres ν_{3}, predicted by the analytical results (see Fig. 1), are visible, especially the resonance 1 : 1, which can lead to chaotic orbits for high enough values of the libration amplitude. The horseshoeshaped orbits at the bottom are mostly chaotic, as expected, since m_{1} + m_{2} = 2 × 10^{−4} is close to the limit 3 × 10^{−4} of their existence (Leleu et al. 2015). The main branch around the 1:1 resonance between ν and ν_{3} is tidally attractive (see Fig. 7). Systems undergoing tides entering this zone can either converge towards the top of the figure or can become completely chaotic. 
3.2 PseudoFixed Points and Linearisation in Their Vicinity
The total differential system that we consider for our model is the one derived from the Hamiltonian (25), F_{0}, to which we now add the tidal perturbations (Eq. (56)). We denote it F : ℝ^{13} ↦ ℝ^{13}. Here we want to find the equilibria of F and to study the linearised dynamics in their vicinity. However, although F_{0} has equilibria (see Table 1), F has none. The five lines of Eq. (56) corresponding to , and cannot all vanish if D_{j} ≠ 0. Since F_{0} has no equilibria at D_{j} = 0 and does not contribute to these five lines, we conclude that F has no equilibria. If the planets are all synchronised, that is, if the all vanish^{8}, then and (57)
We call pseudoequilibrium of F, or pseudofixed point of F, a point X ∈ ℝ^{13} such that . Even though it does not have equilibria, F has pseudoequilibria, and we find them using an extension of the NewtonRaphsonbased algorithm that we developed in Sect. 2.2.3.
Equation (57) shows that, on a branch of pseudo equilibria of F, the parameter γ (and thus the parameter δ, see Eq. (38)) drifts at a speed proportional to the square of the eccentricities. This means that, with tides, the system travels along the main branch from right to left in Fig. 1 (Delisle et al. 2014) much more quickly when δ > 0 than when δ < 0 (due to high eccentricities for positive δ). As the system travels, whether or not it stays close to the main branch or moves away depends on the linear stability of the differential system F in the vicinity of the branch. That is, it depends on the real parts of the eigenvalues of the linear system associated with F. Since γ is not constant at the pseudofixed points, but drifts at a speed given by Eq. (57), computing the eigenvalues of the linearised system makes sense only if γ drifts slowly enough, that is, only if (58)
where the λ_{k} are the eigenvalues of the linearised system. Indeed, is the timescale of evolution of γ, while (max_{k≤13} ℛ_{e}(λ_{k}))^{1} is the timescale of tidal evolution. When the criterion (58) is fulfilled, γ can be considered constant on the timescale of tidal evolution, and the real parts of the linearised system have physical meaning.
Branch 3 always has high eccentricities (see Fig. 2) and exists only for δ > 5.997. The drift in δ towards negative values is quick at high eccentricity (see Eq. (57)), and so branch 3 is tidally very unstable and uninteresting to us. Branch 2 has low values of the eccentricities at large δ, but the existence of a hyperbolic zone at 5.55 ≤ δ ≤ 5.80 (see Table 1) makes it uninteresting as well since the drift ensures that this zone is reached. Hence, we limit the study of tidal dissipation to the main branch (branch 1).
In Fig. 7 we plot the real parts of the eigenvalues of the linearised system associated with F, along its main branch of pseudo equilibria, which is a little perturbation of the main branch of equilibria of Eq. (25). To guarantee that the condition (58) is well respected, we limit ourselves to δ < 1. This is not really a restriction since the tides ensure that this region is quickly reached. We also plot in the same figure the value of for comparison. Only the eigenvalues that are the perturbations of Eqs. (41) and (42) are plotted. In the absence of a third planet, Couturier et al. (2021) show that the eigenvalue responsible for the exponential increase in the libration amplitude of ξ, and thus for the destruction of the coorbital motion, has a real part
and we normalise by this quantity in Fig. 7 and Appendix F. The region –5.64 ≤ δ ≤ –5.47, at the 1 : 1 secular resonance between the libration frequency of the coorbitals, ν, and the frequency of all the pericentres at the pseudo equilibria, ν_{3}, is such that all the eigenvalues of the system have negative real parts, and we expect this region to be linearly stable^{9}. We show in Sect. 4 that this is indeed the case. The linear stability is only temporary though, since the drift in δ ensures that this region is eventually left. We call this the linearly stable region in the remainder of this work. The range in δ corresponding to the linearly stable region strongly depends on m_{3}/ (m_{1} + m_{2}). In Fig. 8 we display its position in the plane (m_{3}, δ).
Fig. 7 Real parts of the eigenvalues of the linearised system associated with F in the vicinity of its main branch of pseudo equilibria, for −7 ≤ δ ≤ 1 (left) and −6.2 ≤ δ ≤ −5 (right). The planetary masses are as in Table 1 and the resonance chain is 1:1:2. Only the four eigenvalues associated with the four degrees of freedom of the conservative system are represented. The other eigenvalues (like those associated with the rotation rates ω_{j}) are of no interest. The real parts behave erratically at the 1:1 secular resonance between ν and ν_{3} (see Fig. 1) in such a way that all of them are negative for −5.64 ≤ δ ≤ −5.47, yielding a linear stability in this region. For δ ≤ −0.104, at least three out of four real parts are negative. In the region δ ≤ 1 we have , and the criterion (58) is very well respected. Similar figures with other planetary masses are available in Appendix F. 
4 Numerical Simulations and Discussions
In this section we investigate the ability of our model to predict the behaviour of a system in the p : p : p + 1 resonance under tidal dissipation. We also check the results drawn in Sect. 3.2 on the linearised dynamics in the vicinity of the main branch.
4.1 Procedure
We numerically integrate two different sets of equations. The first set, our model, is the differential system F, that is, the vector field derived from Eq. (25) to which we add the tidal perturbations (Eq. (56)). The second set is a direct Nbody simulation of the complete system, with the constantΔt model, given by the set of Eq. (E.1).
When the third planet is absent, the relevant parameters to consider to predict the destruction time of a system of two coorbital planets are (Couturier et al. 2021) (60)
namely the total dissipation rate, the mass ratio, the dissipation rate ratio, and the total coorbital mass ratio. Since we are interested in comparing the lifetimes of the coorbitals when they are inside the resonance chain1 : 1 : 2with their lifetimes when they are alone, we make use of these parameters. Due to the complexity of the dynamics of the resonance chain p : p : p + 1, we do not have analytical expressions depending on the parameters of the lifetime of the system, and trying to draw a complete picture would require a very large number of simulations. Instead, we are interested in performing a small number of simulations with parameters that we judge interesting. Thus, we only show the evolution of two systems for the chain 1 : 1 : 2, whose parameters are given in Table 2. Nevertheless, we performed additional simulations with different choices for the planetary masses and the initial δ. The most interesting ones are presented in Appendix F, where we thoroughly discuss the influence of a larger or smaller value for m_{3}.
For the systems listed in Table 2, the chosen value of δ is such that the beginning of the simulation is at the rightmost point of the linearly stable region; these regions are −5.67 ≤ δ ≤ −5.52 for system 0 and −5.64 ≤ δ ≤ −5.47 for system 1 (see Fig. 7). These systems are thus expected to be initially very stable until they leave this region (due to the drift in δ, see Sect. 3.2).
To integrate the two systems with the simplified model F (Eqs. (25) and (56)), we find, for the given value of δ and the planetary masses, the position of the fixed point of the Hamiltonian (25), and use it as initial condition for the integration, with a shift Δξ = 0.1° in ξ, in order to not start exactly at the fixed point. The pseudo fixed point of the model with tides is very close to the fixed point of Eq. (25), and we ignore the difference. To integrate the system with the direct Nbody set of Eq. (56), we find, for the given value of δ and the planetary masses, the position of the libration centre with the algorithm described in Sect. 2.2.4, and use it as the initial condition for the integration, again with the shift Δξ = 0.1°.
When the coorbital planets are alone, the positivity of ℛ_{e} (λ) in Eq. (59) ensures that the system systematically reaches the horseshoeshaped orbits and is destroyed by close encounters. In this case the time τ_{hs} needed to reach the horseshoeshaped orbits is (Couturier et al. 2021) (61)
where T = 2π/η is the coorbital period. Denoting τ_{dest} the coorbital lifetime without third planet, τ_{hs} does not significantly differ from τ_{dest} for a wide range of total coorbital mass, and as long as 10^{−9} ≲ ε ≲ 0.005, we have 1/2 ≲ τ_{dest}/τ_{hs} ≲ 2 (Couturier et al. 2021). We can thus consider that τ_{hs} is the lifetime of the coorbital pair, in the absence of the third planet^{10}. For Δξ = 0.1°, in the case of two coorbital Earthlike planets (using the tidal parameters of Lainey 2016), we have^{11} (62)
We give in Appendix G the time τ_{hs}, computed from Eq. (61), for a variety of hypothetical coorbital pairs. We normalise the time by τ_{hs} in Figs. 9, 10, 11, and in Appendix F.
Parameters of the two numerical simulations.
Fig. 8 Position of the linearly stable region for the resonance chain 1:1:2 in the plane (m_{3}, δ). For every point in this plane the eigenvalues of the linearised system associated with the vector field F are computed, and the point is plotted only if all the real parts are negative. The coorbital masses are m_{1} = m_{2} = 10^{−4} and the tidal parameters are those of the system 0 in Table 2. The position of the linearly stable region weakly depends on m_{1}/m_{2} and on the tidal parameters. The colour gives n_{1}/n_{3} and shows that the chain stabilises the dynamics far from the Keplerian resonance (for which n_{1}/n_{3} = 2). The dashed yellow line plots the secular 1 : 1 resonance between ν and ν_{3}, computed with Eqs. (42) and (33), respectively. For m_{3} > 18 (m_{1} + m_{2}), the linearly stable region disappears, while for m_{3} < 0.29 (m_{1} + m_{2}), two distinct linearly stable regions exist, whose widths tend to 0 with m_{3} (see Appendix F for a discussion of the impact of m_{3} on the dynamics). 
4.2 Mechanisms of CoOrbital Stabilisation
In Figs. 9 and 10 the angles ξ and σ_{j} are plotted as a function of time. The system spends a large amount of time close to the main branch of equilibria, which allows the coorbitals to live notably longer with the presence of the third planet. This can be seen from the destruction occurring at a time t > τ_{hs}. When the system crosses the linearly stable region, the libration amplitude of ξ decreases instead of increasing exponentially, since the real parts of all the eigenvalues of the linearised system associated with F are negative. When the system leaves this region due to the drift in δ and the real part of one eigenvalue becomes positive again, the libration amplitude of ξ is much smaller than it was before entering the linearly stable region. As a result, the system needs more time to reach large libration amplitudes and settle in horseshoeshaped orbits, which delays the coorbital destruction; in other words, crossing the linearly stable region while being sufficiently close to the main branch (so that the linear dynamics dominates) guarantees a coorbital lifetime longer than without the third planet.
In Fig. 9, the negative real parts of all the eigenvalues in the linearly stable region allow the libration amplitude of ξ to reach values as small as 1.8 arcsec at t = 0.57 τ_{hs}. This minimum happens after the linearly stable region is left since the proper mode associated with the newly positive real part (in purple in Fig. 7) has been completely damped by the linearly stable region and some time is needed to pump it noticeably. Similarly, in Fig. 10, 3.5 arcsec of libration amplitude are reached at t = 1.95 τ_{hs}. The early augmentation of the libration amplitude of the angles in the bottom plot is due to the fact that in the direct simulation, the linearly stable region is not exactly at the same values of δ as in the simplified model.
The linearly stable region is not the only reason why the coorbitals in resonant chains can live longer. Another phenomenon, which we refer to as eccentricity damping stabilisation in the rest of this work, allows the libration amplitude to not cross the separatrix leading to horseshoeshaped orbits. After an exponential increase in the libration amplitudes of ξ and σ_{j}, due to at least one eigenvalue with a strictly positive real part, the amplitudes suddenly decrease and the system returns close to the equilibria. This stabilisation of ξ, due to eccentricity damping (see Fig. 11), can happen several times before horseshoeshaped orbits are finally reached, and the system is destroyed (see Fig. 10). The explanation of the eccentricity damping stabilisation relies on the behaviour of the eccentricities. In Fig. 11, we plot the eccentricity e_{1} of planet 1 as a function of time, together with a schema of its behaviour in the plane (e cos ϖ, e sin ϖ). At time t = t_{1} the system is still close to the fixed points, hence the quantity (or any of the two other eccentricities) describes a circle of small radius. Outside the linearly stable region, the eigenvalues of the linearised systems have one positive real part, and as time evolves the radius of the circle grows, while its centre (the equilibrium position of ) gets closer to the origin due to the drift in δ. At time t > t_{2} the circle surrounds the origin and keeps growing, as predicted by the eigenvalues, which triggers a jump in the eccentricity. On the one hand, the linearised system predicts that the circle drawn by grows to infinity; on the other hand, tides impose an exponential decay of the eccentricities. Indeed, the first line of Eq. (56) yields (Correia 2009) (63)
and so, at time t = t_{3}, the circle reaches its maximum radius, the system is now far from its equilibrium, and tides, through nonlinear contributions of the vector field F, force the eccentricities to decrease, which brings the system back to the vicinity of the fixed point at t = t_{4}, where the libration amplitude of the σ_{j}, but also of ξ, is small. In Figs. 9 and 10, several occurrences of the eccentricity damping stabilisation prevent the angle ξ from reaching the horseshoeshaped orbits and increase the lifetime of the coorbital planets.
While the stability induced by the linearly stable region comes from linear contributions of the vector field F, the eccentricity damping stabilisation comes from nonlinear contributions. This latter mechanism works thanks to a strong coupling between the eccentricities (D_{j}, σ_{j}·) and the coorbital angle (L, ξ). In the region δ < 0 (tidally interesting), at most one eigenvalue has a positive real part (see Fig. 7), but due to the coupling, it allows an exponential growth of the libration angle ξ, as well as the eccentricities, which makes the eccentricity damping stabilisation possible. When the time t = t_{3} is reached, the coupling ensures that the eccentricity damping also induces a damping of the libration angle ξ, hence the stabilisation of the coorbital motion. In the absence of the third planet, Couturier et al. (2021) showed that the eccentricities are uncoupled from the coorbital angle ξ. This means that the positive real part (Eq. (59)) associated with (L, ξ) does not induce an exponential growth of the eccentricities, which are on the contrary damped to 0 due to the negative real parts of their eigenvalues. In this case, because of the decoupling, even if some other mechanism increases the eccentricities, the eccentricity damping predicted by Eq. (63) still occurs, but it does not induce a stabilisation of ξ.
The occurrence of the eccentricity damping stabilisation is not systematic. It occurs only if the time t = t_{3} happens before the coorbital planets reach horseshoeshaped orbits. If not, the coorbitals are destroyed before the exponential decrease in the eccentricities can save them. Deciding whether or not a given system will be saved by the eccentricity damping stabilisation requires knowing the proper modes of the linearised system associated with F and how (L, ξ) and the (D_{j}, σ_{j}) are written in the corresponding diagonal basis. Only a numerical work is possible, and we did not undertake it since it is much easier to simply run the corresponding simulation.
If a larger initial δvalue is chosen in these simulations, the system initially has at least one positive real part and moves away from the fixed point at exponential speed. If the eccentricity damping stabilisation works, or if the initial δvalue is small enough, the linearly stable region is reached. However, if the system reaches the linearly stable region when it is too far from the equilibria, the nonlinear contributions of F, combined with the chaotic motion induced by the 1 : 1 secular resonance between ν and ν_{3} (see the stability map in Fig. 6), can lead to peculiar orbits (e.g. switching between the Lagrangian equilibria L_{4} and L_{5}, that is, permutation of the coorbitals). Entering the linearly stable region while still close enough to the equilibria ensures a convergence towards the main branch, and thus an increased stability.
Fig. 9 Evolution of ξ and the σ_{j} (in degrees) as a function of time for system 1 (see Table 2) as integrated by the simplified model F (Eqs. (25) and (56)) (top) and the direct Nbody simulation (Eq. (E.1) (bottom). In the direct simulation horseshoeshaped orbits are reached after 6.42 τ_{hs}, and the coorbital motion is destroyed shortly after that (chaoticity in Fig. 6). The model reaches the horseshoeshaped orbits at t = 13.1 τ_{hs}. Here, the presence of the chain increases the lifetime of the coorbitals by a factor of 6.42. The thickness of the lines in the bottom plot (see e.g. σ_{3}) is due to the shortperiod oscillations that were averaged out in the model. The greyshaded area is the linearly stable region. 
Fig. 10 Evolution of ξ and the aj (in degrees) as a function of time for system 0 as integrated by the simplified model F (Eqs. (25) and (56)) (top) and by the direct Nbody simulation (Eq. (E.1)) (bottom). In the direct simulation the horseshoeshaped orbits are reached after 7.8 τ_{hs}, and the coorbital motion is destroyed shortly after that, (chaoticity in Fig. 6). The model reaches horseshoeshaped orbits at 8.1 τ_{hs}. The greyshaded area is the linearly stable region. 
Fig. 11 Value of e_{1} as a function of time for system 1 in Table 2, as integrated by the simplified model F (Eqs. (25) and (56)) (bottom), and a schematic representation of the eccentricity vector e_{1}e^{iϖ1}, for four particular times (top). The schema explains the different stages of the eccentricity damping stabilisation mechanism. The drift in δ brings the centre of the circle drawn by e_{1}e^{iϖ1} (equilibrium value of e_{1} along the main branch) closer to the origin (see Fig. 1). 
4.3 Discussion
It can be seen in Fig. 5 that for δ = –5.46 the system is already far from the exact resonance. For system 1 at δ = –5.46, we have n_{1} /n_{3} = 2.072. As time goes by, δ drifts towards more negative values, and at t = 6.42τ_{hs}, when it is about to reach horseshoeshaped orbits and be destroyed, system 1 verifies δ = –13.77 and n_{1}/n_{3} = 2.186. Similar considerations are valid for system 0, which means that the system is already outside the resonance, but it is still influenced by the chain. As the system leaves the resonance due to the drift in δ, the coupling between the eccentricities and ξ becomes weaker, meaning that the eccentricity damping stabilisation ends up failing. This prevents the coorbitals from living forever. We nevertheless checked that it can work for simulations starting at n_{1} /n_{3} = 2.37. Only positive values of δ allow for ni/n_{3} a value close to (p + 1) /p = 2 (see Fig. 5). When a positive value of δ is chosen at t = 0, the quick drift in δ due to the high values of the eccentricities (see Eq. (57)) forces the system to reach the region δ < 0 on a timescale much shorter than the timescale of the increase in the libration amplitude of the angles. For systems on the main branch this means that tides favour ratio n_{1} /n_{3} values above their Keplerian value. This result was shown by Delisle et al. (2014) for a twoplanet chain and is confirmed by the observations of the Kepler mission, where a large number of exoplanets were discovered with a mean motion ratio slightly higher than (p + 1) /p (e.g. Delisle & Laskar 2014).
This section shows that our simplified model (Eqs. (25) and (56)) is able to satisfyingly predict the tidal evolution of a resonance chain of the form p : p : p + 1, at least qualitatively, since a precise quantitative description can only be achieved by running the simulation of the direct set of Eq. (E.1). This contrasts with the analytical work performed by Couturier et al. (2021) in the case of lone coorbitals, where the secular model is able to quantitatively predict the outcome of the direct simulations of the complete system with less than 1 % relative error (see their Table 3).
The influence of m_{3}/ (m_{1} + m_{2}) on the coorbital dynamics is thoroughly discussed in Appendix F. We show that, for a large m_{3}, the linearly stable region is not efficient in stabilising the libration amplitude of ξ, while the eccentricity damping stabilisation is very efficient (see Fig. F.1). As m_{3} decreases, the eccentricity damping stabilisation loses efficiency until it does not occur for very small m_{3}values (see Fig. F.3). The linearly stable region has a maximum efficiency for m_{3} ≈ 0.29 (m_{1} + m_{2}) (see Fig. F.2). For a small value of m_{3}, it can stabilise ξ only in a tiny neighbourhood around the 1:1 secular resonance between ν and ν_{3} (see Fig. F.3). Finally, a premature destruction of the coorbital motion (at t < τ_{hs}) can occur for a very large m_{3}value, if the system is already far from the resonance, at a δvalue much lower than its value at the 1 : 1 secular resonance between ν and ν3 (see Fig. F.4).
In terms of coorbital lifetime, the worstcase scenario occurs when the eccentricity damping stabilisation does not work and when the linearly stable region is not crossed (or is crossed while too far from the main branch). In these cases the only positive real part in the region δ < 0 often has a value close to ℛ_{e} (λ) (Fig. 7), and the system reaches the horseshoeshaped orbits at a time close to τ_{hs}.
In summary, in most cases the resonance chain increases the coorbital lifetime, but in a few cases it can also decrease it, especially when m_{3} >> m_{1} + m_{2} and δ is very negative (see Fig. F.4).
5 Conclusion
In this work, we have studied the dynamics of a pair of coorbital planets in the presence of a firstorder resonance with a third planet orbiting outside the coorbitals. We have shown that for systems deep inside this resonance, many equilibria (or libration centres) exist, at least three of which are stable. The existence of a secular resonance between the libration frequency of the coorbitals and the precession frequency of the pericentres can lead to chaotic orbits in the conservative case. However, when tides are involved, we show that this resonance stabilises the coorbital dynamics. Another stabilisation mechanism, due to eccentricity damping, is presented and explained. The model that we built is able to predict the position of the libration centres of the complete system, and we developed an algorithm to find them exactly. When tides are involved, the model reliably gives the qualitative behaviour of the system and, to a certain extent, its quantitative behaviour.
This work shows that when tidal dissipation is included, coorbital systems are more stable if they are inside a resonance chain of the form p : p : p + 1, which increases the chances of a stilltocome detection of a coorbital pair of exoplanets, since Leleu et al. (2019) have shown that coorbital pairs are often formed within a resonance chain. While the analytical work of this paper is performed for any value of the integer p, figures are restricted to the chain 1 : 1 : 2 where p = 1. We nevertheless checked that higher values of p do not impact the qualitative results.
One important contribution of this work is the discovery of a 1 : 1 secular resonance between the libration of the critical angle λ_{1}  λ_{2} and the precession of the pericentres ϖ_{j}, as well as the inherent dynamical consequences (see Figs. 6 and 7). Libration centres, quasiperiodic orbits of the unaveraged problem that generalise the equilibria of the averaged model (Sect. 2.2.1), are such that all the pericentres precess at the same frequency in their vicinity, which holds true for every resonance chain of any number of planets. We thus expect, for other resonance chains, the existence of similar secular resonances, for example between ξ = λ_{1}  4λ_{2} + 3λ_{3} and the mj for the chain 1:2:3.
Acknowledgments
Acknowledgements.
The authors thank JeanBaptiste Delisle and Adrien Leleu for fruitful discussions and subsequent improvements of the paper. This work was supported by CFisUC (UIDB/04564/2020 and UIDP/04564/2020), GRAVITY (PTDC/FISAST/7002/2020), PHOBOS (POCI010145FEDER029932), and ENGAGE SKA (POCI010145FEDER022217), funded by COMPETE 2020 and FCT, Portugal.
Appendix A Notations
For convenience, we list in Table A.1 the notations used throughout this work^{12}.
Notations used in this paper
Appendix B Coefficients in the Expansion of the Hamiltonian
We give here the expressions of the coefficients appearing in Eqs. (21) and (22). They depend on the Laplace coefficients (α) (Laskar & Robutel 1995), and to improve readability we denote (α), where α = ā/a_{3},_{0}. For the resonance 1:1:2 we have (B.1)
for the first order in eccentricity and (B.2)
for the second order.
Appendix C Simplified Differential System
We give in this appendix the expression of the differential system derived from the Hamiltonian ℋ_{K} + ℋ^{(0)} + ℋ^{(1)}which is the Hamiltonian (25) truncated at first order in eccentricity, after the simplifications stated in Sect. 2.2.2
have been performed. We have (C.1)
Appendix D Expression of the matrix Q_{6}
We give in this appendix the matrix Q_{6} appearing in Eq. (40). We denote (D.1)
Appendix E Direct Complete Model for Tides
The complete equations of motion governing the tidal evolution of a planar (N + 1)body system in a heliocentric reference frame, using a linear constant timelag tidal model, are given, for 1 ≤ j ≤ N, by (Mignard 1979) (E.1)
where r_{j} is the heliocentric position vector, θ_{j} the rotation angle of the planet j, k is the unit vector normal to the orbital plane, and f_{j} is the force arising from the tidal potential energy created by the deformation of planet j (Eq. (46)): (E.2)
Appendix F More Numerical Simulations
In this appendix we present the six most interesting simulations that were not shown in Sect. 4. We particularly focus on the influence of the mass m_{3} on the coorbital dynamics. All the simulations comply with m_{1} = m_{2} = 10^{4} m_{0}, and their tidal parameters are those of system 0 in Table 2. We only integrate here the simplified model F (Eq. (56)) (see Sect. 3.2). For each simulation, the real parts of the eigenvalues of the linearised system associated with F are shown alongside the time evolution of the angles ξ and σ_{j}. In the figures of the real parts, a dashed vertical black line shows the starting value of δ, denoted δ_{0}, of the corresponding simulation. In the figures of the angles, a greyshaded area shows the linearly stable region, when relevant. Choosing other tidal parameters does not significantly modify the figures shown here, since we normalise the real parts by ℛ_{e} (λ) (see Eq. (59)) and the times by τ_{hs} (see Eq. (61)), which is the time to reach horseshoeshaped orbits (close to the destruction time) in the absence of a third planet. We invite the reader to have a look at Fig. 8 before reading this appendix, as m_{3} and δ are the only varying parameters between the different simulations. For each δ_{0} and m_{3}, the simulation starts at the corresponding point of the main branch, with a shift Δξ = 0.1° to ξ, in order not to start exactly at equilibrium.
In Fig. F.1, we have m_{3} = 8 (m_{1} + m_{2}), which corresponds to the value yielding the largest linearly stable region (see Fig. 8). However, for this choice of m_{3}, the linearly stable region is located at larger δvalues, and so it drifts quickly (see Sect. 3.2). Furthermore, Eq. (57) shows that the drift in δ is proportional to , and so, due to the rather large m_{3}value, the system leaves the linearly stable region after less than 0.1 τ_{hs}, much more quickly than in Fig. F.2, where m3 is smaller. The amplitude of ξ reaches 0.07° at its lowest, not significantly smaller than its initial value of 0.1°. However, for this choice of m_{3}, the eccentricity damping stabilisation is very efficient, yielding a coorbital lifetime of 8 τ_{hs}.
In Fig. F.2, we have m_{3} = 0.29 (m_{1} + m_{2}), which is the best compromise between the width of the linearly stable region and speed of the drift in δ. For this value of m3, the linearly stable region splits into two distinct strips (see Fig. 8). The system stays in this region for 2 τ_{hs}, 20 times longer than in Fig. F.1, and at t = 3.2 τ_{hs}, ξ librates with only 0.0000012° of amplitude, gaining a factor of 80 000 from the initial 0.1°. The system stays close to the main branch for 6 τ_{hs}, longer than any other simulation that we performed. The eccentricity damping stabilisation is not as efficient as for m_{3} = 8 (m_{1} + m_{2}), but still allows the system to reach 8.5 τ_{hs} before the destruction of the coorbital motion.
Fig. F.1 Evolution of ξ and the σ_{j} for m_{3} = 8 (m_{1} + m_{2}). 
Fig. F.2 Evolution of ξ and the σ_{j} for m_{3} = 0:29(m_{1} + m_{2}). 
In Fig. F.3, we have m_{3} = (m_{1} + m_{2}) /32. With such a low mass for the third planet, the two linearly stable regions are extremely narrow, and the linear stability occurs only if the system is exactly at the 1 : 1 secular resonance between ν and ν_{3} (the initial semimajor axes between the two simulations are very close with a tiny difference in δ). Otherwise, the unique positive real part is ℛ_{e} (λ) everywhere and the destruction occurs at τ_{hs}, as in the absence of the third planet, since the eccentricity damping stabilisation does not work.
Fig. F.3 Evolution of ξ and the σ_{j} for m_{3} = (m_{1} + m_{2}) /32. 
In Fig. F.4, the third planet has a mass m_{3} = 19 (m_{1} + m_{2}) and the linearly stable region does not exist (see Fig. 8). We perform a simulation at the 1:1 secular resonance between ν and ν_{3}, where it would be if it existed. In this region the positive real parts are not greater than ℛ_{e} (λ), and with the eccentricity damping stabilisation, the coorbitals almost reach 4 τ_{hs}. However, very far from the Keplerian resonance (n_{1}/n_{3} = 2.6 at δ = 6, for this choice of m_{3}), the unique positive real part is significantly greater than ℛ_{e} (λ). Furthermore, as (L, ξ) and (D_{j}, σ_{j}) are weakly coupled far from the Keplerian resonance, the eccentricity damping stabilisation is not efficient, leading to the premature destruction (here 0.8 τ_{hs}) of the coorbital motion. Choosing a smaller δ leads to even quicker destruction.
Fig. F.4 Evolution of ξ and the σ_{j} for m_{3} = 19 (m_{1} + m_{2}). 
Appendix G τ_{hs} for hypothetical coorbital pairs
We give in Table G.1 the time Ths for hypothetical coorbital pairs of exoplanets made up of Solar System bodies. The semimajor axis is a = 0.04 AU and the mass of the host star is m_{0} = m_{⊙}, but T_{hs} is easily deduced for other values using the exponents of Eq. (62). We choose Δξ = 0.1°, and again, it is straightforward to extend the results to another choice of Δξ, since τ_{hs} α In (60°/Δξ) (see Eq.(61)).
τ_{hs} for some coorbital systems.
In this table, the tidal parameters are those of Lainey (2016) and only the five bodies for which κ_{2} /Q is well constrained are included. The close Ths between some systems is purely coincidental. It is due to the particular value of Jupiter's κ_{2}/Q and to the fact that this body is much larger and much more massive than the other four.
References
 Beaugé, C., Michtchenko, T.A., & FerrazMello, S. 2006, MNRAS, 365, 1160 [Google Scholar]
 Correia, A. C. M. 2009, ApJ, 704, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Couetdic, J., Laskar, J., Correia, A. C. M., Mayor, M., & Udry, S. 2010, A&A, 519, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Couturier, J., Robutel, P., & Correia, A. C. M. 2021, Celest Mech Dyn Astr, 133, 37 [NASA ADS] [CrossRef] [Google Scholar]
 Cresswell, P., & Nelson, R. P. 2009, A&A, 493, 1141 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Darwin, G. H. 1880, Philos. Trans. R. Soc. Lond. Ser. I, 171, 713 [Google Scholar]
 Delisle, J.B. 2017, A&A, 605, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Delisle, J. B., & Laskar, J. 2014, A&A, 570, L7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Delisle, J.B., Laskar, J., & Correia, A. C. M. 2014, A&A, 566, A137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Deprit, A. 1969, Celest. Mech., 1, 12 [Google Scholar]
 Euler, L. 1764, Novi commentarii academiae scientiarum Petropolitanae. Berlin Acad., 10, 544 [Google Scholar]
 Gascheau, G. 1843, C. R. Acad. Sci. Paris, 16, 393 [Google Scholar]
 Giuppone, C. A., Beaugé, C., Michtchenko, T. A., & FerrazMello, S. 2010, MNRAS, 407, 390 [NASA ADS] [CrossRef] [Google Scholar]
 Henrard, J., & Lemaitre, A. 1983, Celest. Mech., 30, 197 [Google Scholar]
 Kaula, W. M. 1964, Rev. Geophys. Space Phys., 2, 661 [Google Scholar]
 Lagrange. 1772, Œuvres complètes (GouthierVillars, Paris (1869)) [Google Scholar]
 Lainey, V. 2016, Celest. Mech. Dyn. Astron., 126, 145 [Google Scholar]
 Laskar, J. 1990, Icarus, 88, 266 [Google Scholar]
 Laskar, J. 1993, Physica D Nonlinear Phenomena, 67, 257 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J., & Robutel, P. 1995, Celest. Mech. Dyn. Astron., 62, 193 [Google Scholar]
 Laughlin, G., & Chambers, J. E. 2002, Astron. J., 124, 592 [NASA ADS] [CrossRef] [Google Scholar]
 Leleu, A., Robutel, P., & Correia, A. C. M. 2015, A&A, 581, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Leleu, A., Robutel, P., & Correia, A. C. M. 2018, Celest. Mech. Dyn. Astron., 130, 24 [NASA ADS] [CrossRef] [Google Scholar]
 Leleu, A., Coleman, G. A. L., & Ataiee, S. 2019, A&A, 631, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mignard, F. 1979, Moon Planets, 20, 301 [Google Scholar]
 Munk, W. H., & MacDonald, G. J. F. 1960, The Rotation of the Earth; A Geophysical Discussion [Google Scholar]
 Namouni, F. 1999, Icarus, 137, 293 [NASA ADS] [CrossRef] [Google Scholar]
 Pucacco, G. 2021, Celest. Mech. Dyn. Astron., 133, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Robutel, P., & Gabern, F. 2006, MNRAS, 372, 1463 [NASA ADS] [CrossRef] [Google Scholar]
 Robutel, P., & Pousse, A. 2013, Celest. Mech. Dyn. Astron., 117, 17 [Google Scholar]
The exact value of the lower bound is 2 arcsin 23.9° (e.g. Robutel & Pousse 2013).
Even though it is slightly chaotic along the main branch for the conservative system, see the map in Fig. 6.
All Tables
All Figures
Fig.1 Values of the eccentricities and fundamental frequencies given by the analytical results. Left: value of e_{1} at the fixed points of the resonance chain 1:1:2 as predicted by Eq. (35). Right: values of ν/η and ν_{3}/η along the main branch as a function of δ, for the same resonance chain, predicted by Eqs. (42) and (33). In both panels the planetary masses are (m_{1} + m_{2}) /2 = m_{3} = 10^{−4} m_{0} and m_{2}/m_{1} = 10. According to Eq. (32), e_{1} = e_{2}, and for this choice of masses and resonance chain e_{3}/e_{1} = 0.5471. In Table 1, which gathers the fixed points at second order in eccentricity, the three branches visible in the left plot are branches 1 (main branch), 2, and 6. The secular resonances between ν and ν_{3} are shown on the right, and are also easy to spot in Fig. 6. 

In the text 
Fig. 2 Position of the elliptic branches 1, 2, and 3 (see Table 1) of equilibria of Hamiltonian (25) in the resonance 1 : 1 : 2, for −2 ≤ δ ≤ 9. For each branch three curves appear, corresponding to for j = 1,2,3. The planetary masses are the same as in Table 1. A zoomedin image close to the origin is shown in the inset. In this area, the analytical position of these equilibria, given by Eqs. (30) and (35), is plotted by a thin grey line. The agreement is good at low values of eccentricity, but quickly worsens further from the origin. In particular, the thin grey lines are straight since σ_{j}· does not depend on δ in Eq. (30). Branch 1 exists for all values of δ and has all colours from yellow to dark purple, while branch 3 only exists at δ > 5.997 and thus only has purple. 

In the text 
Fig. 3 Trajectories of in the complete Hamiltonian (3) for the six iterations needed for the algorithm to converge to the libration centre. Iteration 0 is the equilibrium (main branch) of the simplified Hamiltonian (25) at δ = 5. It is rather far from the libration centre of the complete Hamiltonian (see also Fig. 4). After six iterations, the algorithm converges to the libration centre and the motion is periodic (hence the closed curves) with frequency ν_{2} (see Sect. 2.2.1). The planetary masses are the same as in Table 1 and the resonance chain is 1:1:2. 

In the text 
Fig. 4 Position of branch 1 of elliptic libration centres of the complete Hamiltonian (3) in the resonance 1:1:2 for −7 ≤ δ ≤ 7. The planetary masses are the same as in Table 1. As a comparison, branch 1 of equilibria of the simplified Hamiltonian (25) is plotted in grey for the same range in δ. For small enough values of δ the eccentricity is not too high and the agreement is good. The simplified and the complete Hamiltonian diverge when δ → +∞. 

In the text 
Fig. 5 ν_{3}/η as a function of δ, in the complete Hamiltonian, for the resonance chain 1 : 1 : 2, along branch 1 of the libration centres in Fig. 4. The analytical expression (33) is plotted in grey for comparison. The colour gives the value of n_{1} /n_{3}. As expected from Eq. (33), ν_{3} is proportional to δ far from the resonance where the eccentricities are small. 

In the text 
Fig. 6 Diffusion index ζ_{ν} as a function of δ and ξ_{0}. The planetary masses are as in Table 1 and the chain is 1:1:2. The top of the figure, at ξ_{0} ≈ 60°, is the main branch. The bluetogreen regions are almost quasiperiodic (stable), while the red regions are chaotic (unstable). The secular resonances between the libration frequency ν and the precession frequency of the pericentres ν_{3}, predicted by the analytical results (see Fig. 1), are visible, especially the resonance 1 : 1, which can lead to chaotic orbits for high enough values of the libration amplitude. The horseshoeshaped orbits at the bottom are mostly chaotic, as expected, since m_{1} + m_{2} = 2 × 10^{−4} is close to the limit 3 × 10^{−4} of their existence (Leleu et al. 2015). The main branch around the 1:1 resonance between ν and ν_{3} is tidally attractive (see Fig. 7). Systems undergoing tides entering this zone can either converge towards the top of the figure or can become completely chaotic. 

In the text 
Fig. 7 Real parts of the eigenvalues of the linearised system associated with F in the vicinity of its main branch of pseudo equilibria, for −7 ≤ δ ≤ 1 (left) and −6.2 ≤ δ ≤ −5 (right). The planetary masses are as in Table 1 and the resonance chain is 1:1:2. Only the four eigenvalues associated with the four degrees of freedom of the conservative system are represented. The other eigenvalues (like those associated with the rotation rates ω_{j}) are of no interest. The real parts behave erratically at the 1:1 secular resonance between ν and ν_{3} (see Fig. 1) in such a way that all of them are negative for −5.64 ≤ δ ≤ −5.47, yielding a linear stability in this region. For δ ≤ −0.104, at least three out of four real parts are negative. In the region δ ≤ 1 we have , and the criterion (58) is very well respected. Similar figures with other planetary masses are available in Appendix F. 

In the text 
Fig. 8 Position of the linearly stable region for the resonance chain 1:1:2 in the plane (m_{3}, δ). For every point in this plane the eigenvalues of the linearised system associated with the vector field F are computed, and the point is plotted only if all the real parts are negative. The coorbital masses are m_{1} = m_{2} = 10^{−4} and the tidal parameters are those of the system 0 in Table 2. The position of the linearly stable region weakly depends on m_{1}/m_{2} and on the tidal parameters. The colour gives n_{1}/n_{3} and shows that the chain stabilises the dynamics far from the Keplerian resonance (for which n_{1}/n_{3} = 2). The dashed yellow line plots the secular 1 : 1 resonance between ν and ν_{3}, computed with Eqs. (42) and (33), respectively. For m_{3} > 18 (m_{1} + m_{2}), the linearly stable region disappears, while for m_{3} < 0.29 (m_{1} + m_{2}), two distinct linearly stable regions exist, whose widths tend to 0 with m_{3} (see Appendix F for a discussion of the impact of m_{3} on the dynamics). 

In the text 
Fig. 9 Evolution of ξ and the σ_{j} (in degrees) as a function of time for system 1 (see Table 2) as integrated by the simplified model F (Eqs. (25) and (56)) (top) and the direct Nbody simulation (Eq. (E.1) (bottom). In the direct simulation horseshoeshaped orbits are reached after 6.42 τ_{hs}, and the coorbital motion is destroyed shortly after that (chaoticity in Fig. 6). The model reaches the horseshoeshaped orbits at t = 13.1 τ_{hs}. Here, the presence of the chain increases the lifetime of the coorbitals by a factor of 6.42. The thickness of the lines in the bottom plot (see e.g. σ_{3}) is due to the shortperiod oscillations that were averaged out in the model. The greyshaded area is the linearly stable region. 

In the text 
Fig. 10 Evolution of ξ and the aj (in degrees) as a function of time for system 0 as integrated by the simplified model F (Eqs. (25) and (56)) (top) and by the direct Nbody simulation (Eq. (E.1)) (bottom). In the direct simulation the horseshoeshaped orbits are reached after 7.8 τ_{hs}, and the coorbital motion is destroyed shortly after that, (chaoticity in Fig. 6). The model reaches horseshoeshaped orbits at 8.1 τ_{hs}. The greyshaded area is the linearly stable region. 

In the text 
Fig. 11 Value of e_{1} as a function of time for system 1 in Table 2, as integrated by the simplified model F (Eqs. (25) and (56)) (bottom), and a schematic representation of the eccentricity vector e_{1}e^{iϖ1}, for four particular times (top). The schema explains the different stages of the eccentricity damping stabilisation mechanism. The drift in δ brings the centre of the circle drawn by e_{1}e^{iϖ1} (equilibrium value of e_{1} along the main branch) closer to the origin (see Fig. 1). 

In the text 
Fig. F.1 Evolution of ξ and the σ_{j} for m_{3} = 8 (m_{1} + m_{2}). 

In the text 
Fig. F.2 Evolution of ξ and the σ_{j} for m_{3} = 0:29(m_{1} + m_{2}). 

In the text 
Fig. F.3 Evolution of ξ and the σ_{j} for m_{3} = (m_{1} + m_{2}) /32. 

In the text 
Fig. F.4 Evolution of ξ and the σ_{j} for m_{3} = 19 (m_{1} + m_{2}). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.