© J. Couturier et al. 2022

1 Introduction

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

In the planar and circular case, for masses in the range 3 × 10⁻⁴ < (m₁ + m₂) /m₀ < 1/27, the angle λ₁ – λ₂ between the co-orbitals librates around its equilibrium of ±60° in trajectories commonly called tadpole orbits, where the co-orbital angle is bounded by¹ 23.9° < λ₁ – λ₂ < 180°. However, for masses (m₁ + m₂) /m₀ < 3 × 10⁻⁴ (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 λ₁ – λ₂ 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 anti-Lagrange have been shown, numerically by Giuppone et al. (2010) and analytically by Robutel & Pousse (2013). In the Lagrange configuration, the pericentres of the orbits do not precess and verify the relation ϖ₁ – ϖ₂ = 60°, whereas in the anti-Lagrange configuration, both orbits precess at the same frequency while maintaining the relation ϖ₁ – ϖ₂ = 240^°. For low eccentricities but very large libration amplitudes, in the region of horseshoe-shaped orbits, Couturier et al. (2021) showed that the Lagrange and anti-Lagrange configurations correspond to ϖ₁ – ϖ₂= 0° and ϖ₁ – ϖ₂ = 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 co-orbitals (Namouni 1999).

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

For co-orbital exoplanets orbiting close to the host star, tidal dissipation induced by the differential gravitational interaction leads to a long-term evolution of the orbits. Couturier et al. (2021) shows that, for a pair of co-orbital exoplanets orbiting a star, the equilateral Lagrangian equilibria are always repulsive under tidal interactions, and that regardless of the parameters of the system, the destruction of the co-orbital motion is unavoidable. However, the discovery of co-orbital exoplanets is still possible because the destruction time is strongly dependent on the parameters of the system and can easily be greater than the lifetime of the host star. Couturier et al. (2021) neglected any interaction with possible other planets in the system. In this paper we extend this work to the case where the pair of co-orbital exoplanets interacts with an outermost third planet, in a first-order mean motion resonance with the co-orbitals. More precisely, we study the Laplace resonance chain p : p : p + 1, where p is a small integer.

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

2 . The p : p : p +1 Resonance Chain

2.1 The Hamiltonian of the Problem

In this section, we study an occurrence of the point-mass planar four-body problem. We construct the Hamiltonian associated with the resonance chain p : p : p + 1, where a central body, a star of mass m₀, is orbited by two co-orbital planets of mass m₁ and m₂ and a third planet of mass m₃ further away from the star and in a first-order mean motion resonance with the pair of co-orbital planets. Although we write all equations for a general value of the integer p, the figures are restricted to the case p = 1 where the nominal period of the third planet is twice that of the co-orbitals. For all planets we define the quantities β_j = m₀m_j/(m₀ + 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 semi-major 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 star-planet interactions, reads (4)

and the perturbation H_P, whose size relative to H_K is of order (m₁ + m₂ + m₃) /m₀, 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 semi-major axes a_j,0 are related to n_j,0 by the Kepler law . The a_j are always close to their nominal value² 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 ξ₃. Moreover, in the p : p : p + 1 resonance, the angle ξ₂ 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 (ξ₂, ξ₃, Γ′, 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₁ + D₂ + D₃ = ∑_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³. 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₁ + D₂ + D₃ = ∑_j ΔΛ_j. We find the Hamil-tonian (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₁, p₂) ϵ {(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 non-zero , the conservation of the angular momentum imposes on the tuples q = (q₁, q₂, q₃, q₄) ϵ ℕ⁴ and k = (k₁, k₂) ϵ ℤ² 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 ℋ⁽⁰⁾ has no subscript since only the pair of co-orbitals yields terms of order 0 and no confusion is possible.

Following Laskar & Robutel (1995), ℋ_1,3and ℋ_2,3can be written (21)

and (22)

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,2cannot be obtained using the same procedure since the Laplace coefficients diverge for equal nominal semi-major axes and a_1,0 = a_2,0 = ā. We instead follow the method described in Robutel & Pousse (2013). We denote and find (23)

where (24)

The final simplified Hamiltonian is then (25)

We denote F₀ : ℝ⁸ ↦ ℝ⁸ the differential system derived from Eq. (25) by the Hamilton-Jacobi equations.

2.2 Equilibria and Linearisation in Their Vicinity

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

2.2.1 Fixed Points and Libration Centres

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

At a fixed point (or equilibrium) of the model, L and D_j are constant and so are e_j and a_j. Similarly, the angles σ_j and ξ are constant; that is, there are constants c_j such that (26)

However, the secular angle –pλ₂ + (p + 1) λ₃ and the pericentres a j are not constant at the equilibria, but they all precess with the same frequency, which we denote ν₃.

The average performed in Eq. (10) is actually analogous to a first-order Lie serie expansion, and the averaging process can be seen as a periodic change of variable. Denoting x and x′ the variables of the Hamiltonian respectively before and after the average, we have , 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λ₂ + (p + 1) λ₃ or the pericentres ϖ_j, a fixed point in the model corresponds in the complete Hamiltonian to a quasi-periodic motion with the two frequencies ν₂ and ν₃, with (28)

where these qua tities are evaluated at the equilibrium. More precisely, it is a periodic motio with frequecy ν₂ i a rotatig frame followig all the perice tres at frequecy ν₃. 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 first-order resonances between non-consecutive planets. However, we show here that a further simplification allows us to obtain analytical expressions of the equilibria and of the eigenvalues of the linearised system.

To further simplify the Hamiltonian, we force a decoupling between the degree of freedom (ξ, ΔL) associated with the libration of the co-orbitals and the three other degrees of freedom (σ_j, D_j) To this end, we first evaluate at ξ = π/3. Indeed, the Hamiltonian ℋ_K + ℋ⁽⁰⁾ + ℋ⁽¹⁾ is only a perturbation of ℋ_K + ℋ⁽⁰⁾ 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 anti-diagonal term ⁴ 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 + ℋ⁽⁰⁾ + ℋ⁽¹⁾ once these simplifications have been performed, is given in Appendix C. It vanishes when the angles are equal to (30)

where (31)

and when the actions are equal to⁵ (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, ν₃, 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 third-degree polynomial in Z:

Z³ – 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² – 4P³ = 0, that is at (36)

In the rest of this work, we use the parameter (37)

where (38)

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 third-degree 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₁,u₂,u₃,υ₁,υ₂,υ₃, ΔL, ξ), the linearised systemis (40)

and X₀ is the equilibrium value of X. The matrix Q₆ 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₃, 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₂ 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₃ in the vicinity of the main branch shows that, at δ ≈ −5.6 for the planetary masses in Fig. 1, all these eigenvalues have roughly the same value, yielding a 1 : 1 secular resonance between the libration frequency of the co-orbitals and the precession frequency of the pericentres. Other secular resonances between ν and ν3 are shown in Fig. 1, and are also very visible on the stability map from Fig. 6. We show in Sect. 3.2 that the secular resonance 1 : 1 has important consequences for the tidal stability of the co-orbital pair.

Fig.1 Values of the eccentricities and fundamental frequencies given by the analytical results. Left: value of e1 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 (m1 + m2) /2 = m3 = 10−4 m0 and m2/m1 = 10. According to Eq. (32), e1 = e2, and for this choice of masses and resonance chain e3/e1 = 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.

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⁶, 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₀ derived from the Hamiltonian (25) depends on the choice of the parameter δ (through ΔG), and once a fixed point is found for a particular value of δ, we repeat the Newton-Raphson algorithm for slowly varying values of δ in order to travel along the whole branch. We look for equilibria exploring the parallelepiped in the phase space defined by|u_j| < 0.08 and . We choose as initial condition of the Newton-Raphson method for L since all equilibria are expected to be close to this value, and so no discretisation is necessary along this axis. In the same way, we only choose ξ₀ ϵ {±π/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⁶ 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 two-planet 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 zoomed-in image close to the origin is shown in the inset. In this area, the analytical position of these equilibria, given by Eqs. (30) and (35), is plotted by a thin grey line. The agreement is good at low values of eccentricity, but quickly worsens further from the origin. In particular, the thin grey lines are straight since σ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 quasi-periodic, and we write, for any complex quantity z depending on these trajectories, (43)

where the coordinates of ω = ^t(ν, ν₂, ν₃, g₁, g₂, g₃) are the fundamental frequencies of the complete Hamiltonian (3). The frequencies ν, ν₂, and ν₃ 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⁷ with frequency ν₂. 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 ν₂ 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 ν₂ 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 ν₂ is much larger than the other frequencies, hence it is easy to isolate terms that depend on ν₂ 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 ν₂.

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, ν₃, 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₁/n₃ 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 ν₃ 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 point-mass p : p : p + 1 resonance chain by constructing a dynamical map using the frequency analysis method to determine the chaoticity of a given orbit (Laskar 1990). More precisely, we study the stability of the chain along its main branch (see Fig. 1) between δ = 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 horseshoe-shaped orbits (close to 24°, see Robutel & Pousse 2013). For values of ξ₀ close to 60°, the considered orbit is close to the main branch and it moves away for decreasing values of ξ₀. For all other variables, we choose as the initial condition the value at the libration centre. Every trajectory (i.e. every choice of δ and ξ₀) 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 ν⁽²⁾ is extracted from the frequency analysis of e^iξ. We obtain two values of ν, namely ν⁽¹⁾ for the first half of 40 000 periods and ν⁽²⁾ for the second half. The diffusion index, defined as (Robutel & Gabern 2006) (44)

measures the degree of quasi-periodicity of the orbit. Orbits with ζ_ν < −6 are considered close to quasi-periodic (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 ν₃, 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 n1 /n3. 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₂ 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)² = 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)

with (51)

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)

and (55)

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 Hamilton-Jacobi equations and considering the starred variables as parameters. The starred variables are then expressed by a first-order Taylor expansion in Eq. (47) (see Couturier et al. 2021 for more details). The perturbation to the vector field, due to tides and at second order in eccentricity, reads (56)

where we posed Since the differential system does not depend on ξ₂ and ξ₃ and as their dynamics are of no interest to us, the lines ξ₂ and ξ₃ 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 blue-to-green regions are almost quasi-periodic (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 horseshoe-shaped orbits at the bottom are mostly chaotic, as expected, since m1 + m2 = 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 Pseudo-Fixed Points and Linearisation in Their Vicinity

The total differential system that we consider for our model is the one derived from the Hamiltonian (25), F₀, to which we now add the tidal perturbations (Eq. (56)). We denote it F : ℝ¹³ ↦ ℝ¹³. Here we want to find the equilibria of F and to study the linearised dynamics in their vicinity. However, although F₀ has equilibria (see Table 1), F has none. The five lines of Eq. (56) corresponding to , and cannot all vanish if D_j ≠ 0. Since F₀ 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⁸, then and (57)

We call pseudo-equilibrium of F, or pseudo-fixed point of F, a point X ∈ ℝ¹³ such that . Even though it does not have equilibria, F has pseudo-equilibria, and we find them using an extension of the Newton-Raphson-based algorithm that we developed in Sect. 2.2.3.

Equation (57) shows that, on a branch of pseudo equilibria of F, the parameter γ (and thus the parameter δ, see Eq. (38)) drifts at a speed proportional to the square of the eccentricities. This means that, with tides, the system travels along the main branch from right to left in Fig. 1 (Delisle et al. 2014) much more quickly when δ > 0 than when δ < 0 (due to high eccentricities for positive δ). As the system travels, whether or not it stays close to the main branch or moves away depends on the linear stability of the differential system F in the vicinity of the branch. That is, it depends on the real parts of the eigenvalues of the linear system associated with F. Since γ is not constant at the pseudo-fixed points, but drifts at a speed given by Eq. (57), computing the eigenvalues of the linearised system makes sense only if γ drifts slowly enough, that is, only if (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 co-orbital motion, has a real part

(59)

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 co-orbitals, ν, and the frequency of all the pericentres at the pseudo equilibria, ν₃, is such that all the eigenvalues of the system have negative real parts, and we expect this region to be linearly stable⁹. 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₃/ (m₁ + m₂). In Fig. 8 we display its position in the plane (m₃, δ).

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 N-body 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 co-orbital mass ratio. Since we are interested in comparing the lifetimes of the co-orbitals 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₃.

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 N-body 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 co-orbital planets are alone, the positivity of ℛ_e (λ) in Eq. (59) ensures that the system systematically reaches the horseshoe-shaped orbits and is destroyed by close encounters. In this case the time τ_hs needed to reach the horseshoe-shaped orbits is (Couturier et al. 2021) (61)

where T = 2π/η is the co-orbital period. Denoting τ_dest the co-orbital lifetime without third planet, τ_hs does not significantly differ from τ_dest for a wide range of total co-orbital mass, and as long as 10⁻⁹ ≲ ε ≲ 0.005, we have 1/2 ≲ τ_dest/τ_hs ≲ 2 (Couturier et al. 2021). We can thus consider that τ_hs is the lifetime of the co-orbital pair, in the absence of the third planet¹⁰. For Δξ = 0.1°, in the case of two co-orbital Earth-like planets (using the tidal parameters of Lainey 2016), we have¹¹ (62)

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

Fig. 8 Position of the linearly stable region for the resonance chain 1:1:2 in the plane (m3, δ). For every point in this plane the eigenvalues of the linearised system associated with the vector field F are computed, and the point is plotted only if all the real parts are negative. The co-orbital masses are m1 = m2 = 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 m1/m2 and on the tidal parameters. The colour gives n1/n3 and shows that the chain stabilises the dynamics far from the Keplerian resonance (for which n1/n3 = 2). The dashed yellow line plots the secular 1 : 1 resonance between ν and ν3, computed with Eqs. (42) and (33), respectively. For m3 > 18 (m1 + m2), the linearly stable region disappears, while for m3 < 0.29 (m1 + m2), two distinct linearly stable regions exist, whose widths tend to 0 with m3 (see Appendix F for a discussion of the impact of m3 on the dynamics).

4.2 Mechanisms of Co-Orbital 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 co-orbitals to live notably longer with the presence of the third planet. This can be seen from the destruction occurring at a time t > τ_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 horseshoe-shaped orbits, which delays the co-orbital destruction; in other words, crossing the linearly stable region while being sufficiently close to the main branch (so that the linear dynamics dominates) guarantees a co-orbital lifetime longer than without the third planet.

In Fig. 9, the negative real parts of all the eigenvalues in the linearly stable region allow the libration amplitude of ξ 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 co-orbitals in resonant chains can live longer. Another phenomenon, which we refer to as eccentricity damping stabilisation in the rest of this work, allows the libration amplitude to not cross the separatrix leading to horseshoe-shaped orbits. After an exponential increase in the libration amplitudes of ξ 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 horseshoe-shaped orbits are finally reached, and the system is destroyed (see Fig. 10). The explanation of the eccentricity damping stabilisation relies on the behaviour of the eccentricities. In Fig. 11, we plot the eccentricity e₁ 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₁ 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₂ 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₃, the circle reaches its maximum radius, the system is now far from its equilibrium, and tides, through non-linear contributions of the vector field F, force the eccentricities to decrease, which brings the system back to the vicinity of the fixed point at t = t₄, 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 horseshoe-shaped orbits and increase the lifetime of the co-orbital planets.

While the stability induced by the linearly stable region comes from linear contributions of the vector field F, the eccentricity damping stabilisation comes from non-linear contributions. This latter mechanism works thanks to a strong coupling between the eccentricities (D_j, σ_j·) and the co-orbital 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₃ is reached, the coupling ensures that the eccentricity damping also induces a damping of the libration angle ξ, hence the stabilisation of the co-orbital motion. In the absence of the third planet, Couturier et al. (2021) showed that the eccentricities are uncoupled from the co-orbital angle ξ. 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₃ happens before the co-orbital planets reach horseshoe-shaped orbits. If not, the co-orbitals are destroyed before the exponential decrease in the eccentricities can save them. Deciding whether or not a given system will be saved by the eccentricity damping stabilisation requires knowing the proper modes of the linearised system associated with F and how (L, ξ) 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 non-linear contributions of F, combined with the chaotic motion induced by the 1 : 1 secular resonance between ν and ν₃ (see the stability map in Fig. 6), can lead to peculiar orbits (e.g. switching between the Lagrangian equilibria L₄ and L₅, that is, permutation of the co-orbitals). Entering the linearly stable region while still close enough to the equilibria ensures a convergence towards the main branch, and thus an increased stability.

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 N-body simulation (Eq. (E.1) (bottom). In the direct simulation horseshoe-shaped orbits are reached after 6.42 τhs, and the co-orbital motion is destroyed shortly after that (chaoticity in Fig. 6). The model reaches the horseshoe-shaped orbits at t = 13.1 τhs. Here, the presence of the chain increases the lifetime of the co-orbitals by a factor of 6.42. The thickness of the lines in the bottom plot (see e.g. σ3) is due to the short-period oscillations that were averaged out in the model. The grey-shaded 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 N-body simulation (Eq. (E.1)) (bottom). In the direct simulation the horseshoe-shaped orbits are reached after 7.8 τhs, and the co-orbital motion is destroyed shortly after that, (chaoticity in Fig. 6). The model reaches horseshoe-shaped orbits at 8.1 τhs. The grey-shaded area is the linearly stable region.

Fig. 11 Value of e1 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 e1eiϖ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 e1eiϖ1 (equilibrium value of e1 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₁ /n₃ = 2.072. As time goes by, δ drifts towards more negative values, and at t = 6.42τ_hs, when it is about to reach horseshoe-shaped orbits and be destroyed, system 1 verifies δ = –13.77 and n₁/n₃ = 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 co-orbitals from living forever. We nevertheless checked that it can work for simulations starting at n₁ /n₃ = 2.37. Only positive values of δ allow for ni/n₃ 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₁ /n₃ values above their Keplerian value. This result was shown by Delisle et al. (2014) for a two-planet chain and is confirmed by the observations of the Kepler mission, where a large number of exoplanets were discovered with a mean motion ratio slightly higher than (p + 1) /p (e.g. Delisle & Laskar 2014).

This section shows that our simplified model (Eqs. (25) and (56)) is able to satisfyingly predict the tidal evolution of a resonance chain of the form p : p : p + 1, at least qualitatively, since a precise quantitative description can only be achieved by running the simulation of the direct set of Eq. (E.1). This contrasts with the analytical work performed by Couturier et al. (2021) in the case of lone co-orbitals, where the secular model is able to quantitatively predict the outcome of the direct simulations of the complete system with less than 1 % relative error (see their Table 3).

The influence of m₃/ (m₁ + m₂) on the co-orbital dynamics is thoroughly discussed in Appendix F. We show that, for a large m₃, 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₃ decreases, the eccentricity damping stabilisation loses efficiency until it does not occur for very small m₃-values (see Fig. F.3). The linearly stable region has a maximum efficiency for m₃ ≈ 0.29 (m₁ + m₂) (see Fig. F.2). For a small value of m₃, it can stabilise ξ only in a tiny neighbourhood around the 1:1 secular resonance between ν and ν₃ (see Fig. F.3). Finally, a premature destruction of the co-orbital motion (at t < τ_hs) can occur for a very large m₃-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 co-orbital lifetime, the worst-case scenario occurs when the eccentricity damping stabilisation does not work and when the linearly stable region is not crossed (or is crossed while too far from the main branch). In these cases the only positive real part in the region δ < 0 often has a value close to ℛ_e (λ) (Fig. 7), and the system reaches the horseshoe-shaped orbits at a time close to τ_hs.

In summary, in most cases the resonance chain increases the co-orbital lifetime, but in a few cases it can also decrease it, especially when m₃ >> m₁ + m₂ and δ is very negative (see Fig. F.4).

5 Conclusion

In this work, we have studied the dynamics of a pair of co-orbital planets in the presence of a first-order resonance with a third planet orbiting outside the co-orbitals. We have shown that for systems deep inside this resonance, many equilibria (or libration centres) exist, at least three of which are stable. The existence of a secular resonance between the libration frequency of the co-orbitals and the precession frequency of the pericentres can lead to chaotic orbits in the conservative case. However, when tides are involved, we show that this resonance stabilises the co-orbital dynamics. Another stabilisation mechanism, due to eccentricity damping, is presented and explained. The model that we built is able to predict the position of the libration centres of the complete system, and we developed an algorithm to find them exactly. When tides are involved, the model reliably gives the qualitative behaviour of the system and, to a certain extent, its quantitative behaviour.

This work shows that when tidal dissipation is included, co-orbital systems are more stable if they are inside a resonance chain of the form p : p : p + 1, which increases the chances of a still-to-come detection of a co-orbital pair of exoplanets, since Leleu et al. (2019) have shown that co-orbital pairs are often formed within a resonance chain. While the analytical work of this paper is performed for any value of the integer p, figures are restricted to the chain 1 : 1 : 2 where p = 1. We nevertheless checked that higher values of p do not impact the qualitative results.

One important contribution of this work is the discovery of a 1 : 1 secular resonance between the libration of the critical angle λ₁ - λ₂ and the precession of the pericentres ϖ_j, as well as the inherent dynamical consequences (see Figs. 6 and 7). Libration centres, quasi-periodic orbits of the unaveraged problem that generalise the equilibria of the averaged model (Sect. 2.2.1), are such that all the pericentres precess at the same frequency in their vicinity, which holds true for every resonance chain of any number of planets. We thus expect, for other resonance chains, the existence of similar secular resonances, for example between ξ = λ₁ - 4λ₂ + 3λ₃ and the mj for the chain 1:2:3.

Acknowledgments

Acknowledgements.

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

Appendix A Notations

For convenience, we list in Table A.1 the notations used throughout this work¹².

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₃,₀. 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 + ℋ⁽⁰⁾ + ℋ⁽¹⁾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)

where (C.2)

Appendix D Expression of the matrix Q₆

We give in this appendix the matrix Q₆ appearing in Eq. (40). We denote (D.1)

and obtain Q₆ =

Appendix E Direct Complete Model for Tides

The complete equations of motion governing the tidal evolution of a planar (N + 1)-body system in a heliocentric reference frame, using a linear constant time-lag tidal model, are given, for 1 ≤ 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₃ on the co-orbital dynamics. All the simulations comply with m₁ = m₂ = 10^-4 m₀, 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 δ₀, of the corresponding simulation. In the figures of the angles, a grey-shaded area shows the linearly stable region, when relevant. Choosing other tidal parameters does not significantly modify the figures shown here, since we normalise the real parts by ℛ_e (λ) (see Eq. (59)) and the times by τ_hs (see Eq. (61)), which is the time to reach horseshoe-shaped orbits (close to the destruction time) in the absence of a third planet. We invite the reader to have a look at Fig. 8 before reading this appendix, as m₃ and δ are the only varying parameters between the different simulations. For each δ₀ and m₃, 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₃ = 8 (m₁ + m₂), which corresponds to the value yielding the largest linearly stable region (see Fig. 8). However, for this choice of m₃, 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₃-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₃, the eccentricity damping stabilisation is very efficient, yielding a co-orbital lifetime of 8 τ_hs.

In Fig. F.2, we have m₃ = 0.29 (m₁ + m₂), 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₃ = 8 (m₁ + m₂), but still allows the system to reach 8.5 τ_hs before the destruction of the co-orbital motion.

Fig. F.1 Evolution of ξ and the σj for m3 = 8 (m1 + m2).

Fig. F.2 Evolution of ξ and the σj for m3 = 0:29(m1 + m2).

In Fig. F.3, we have m₃ = (m₁ + m₂) /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 ν₃ (the initial semi-major 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 m3 = (m1 + m2) /32.

In Fig. F.4, the third planet has a mass m₃ = 19 (m₁ + m₂) and the linearly stable region does not exist (see Fig. 8). We perform a simulation at the 1:1 secular resonance between ν and ν₃, 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 co-orbitals almost reach 4 τ_hs. However, very far from the Keplerian resonance (n₁/n₃ = 2.6 at δ = -6, for this choice of m₃), 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 co-orbital motion. Choosing a smaller δ leads to even quicker destruction.

Fig. F.4 Evolution of ξ and the σj for m3 = 19 (m1 + m2).

Appendix G τ_hs for hypothetical co-orbital pairs

We give in Table G.1 the time Ths for hypothetical co-orbital 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₀ = 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)).

In this table, the tidal parameters are those of Lainey (2016) and only the five bodies for which κ₂ /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 κ₂/Q and to the fact that this body is much larger and much more massive than the other four.

References

All Tables

All Figures

Fig.1 Values of the eccentricities and fundamental frequencies given by the analytical results. Left: value of e1 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 (m1 + m2) /2 = m3 = 10−4 m0 and m2/m1 = 10. According to Eq. (32), e1 = e2, and for this choice of masses and resonance chain e3/e1 = 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 zoomed-in image close to the origin is shown in the inset. In this area, the analytical position of these equilibria, given by Eqs. (30) and (35), is plotted by a thin grey line. The agreement is good at low values of eccentricity, but quickly worsens further from the origin. In particular, the thin grey lines are straight since σ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 n1 /n3. 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 blue-to-green regions are almost quasi-periodic (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 horseshoe-shaped orbits at the bottom are mostly chaotic, as expected, since m1 + m2 = 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 (m3, δ). For every point in this plane the eigenvalues of the linearised system associated with the vector field F are computed, and the point is plotted only if all the real parts are negative. The co-orbital masses are m1 = m2 = 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 m1/m2 and on the tidal parameters. The colour gives n1/n3 and shows that the chain stabilises the dynamics far from the Keplerian resonance (for which n1/n3 = 2). The dashed yellow line plots the secular 1 : 1 resonance between ν and ν3, computed with Eqs. (42) and (33), respectively. For m3 > 18 (m1 + m2), the linearly stable region disappears, while for m3 < 0.29 (m1 + m2), two distinct linearly stable regions exist, whose widths tend to 0 with m3 (see Appendix F for a discussion of the impact of m3 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 N-body simulation (Eq. (E.1) (bottom). In the direct simulation horseshoe-shaped orbits are reached after 6.42 τhs, and the co-orbital motion is destroyed shortly after that (chaoticity in Fig. 6). The model reaches the horseshoe-shaped orbits at t = 13.1 τhs. Here, the presence of the chain increases the lifetime of the co-orbitals by a factor of 6.42. The thickness of the lines in the bottom plot (see e.g. σ3) is due to the short-period oscillations that were averaged out in the model. The grey-shaded 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 N-body simulation (Eq. (E.1)) (bottom). In the direct simulation the horseshoe-shaped orbits are reached after 7.8 τhs, and the co-orbital motion is destroyed shortly after that, (chaoticity in Fig. 6). The model reaches horseshoe-shaped orbits at 8.1 τhs. The grey-shaded area is the linearly stable region.

In the text

Fig. 11 Value of e1 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 e1eiϖ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 e1eiϖ1 (equilibrium value of e1 along the main branch) closer to the origin (see Fig. 1).

In the text

	Fig. F.1 Evolution of ξ and the σj for m3 = 8 (m1 + m2).
In the text

	Fig. F.2 Evolution of ξ and the σj for m3 = 0:29(m1 + m2).
In the text

	Fig. F.3 Evolution of ξ and the σj for m3 = (m1 + m2) /32.
In the text

	Fig. F.4 Evolution of ξ and the σj for m3 = 19 (m1 + m2).
In the text