© ESO 2019

1. Introduction

Among the known multiplanetary systems, a significant number contain bodies in (or close to) first and second order mean-motion resonances (MMR) (Fabrycky et al. 2014). However, no planets were found in a zeroth order MMR, also called trojan or co-orbital configuration, despite several dedicated studies (Madhusudhan & Winn 2009; Janson 2013; Lillo-Box et al. 2018a,b). The formation of the first and second order resonances is generally explained by the convergent migration of two planets under the dissipative forces applied by the protoplanetary disc (see for example Lee & Peale 2002). In the co-orbital case, the resonance is surrounded by a chaotic area due to the overlap of first-order MMRs (Wisdom 1980; Deck et al. 2013). The crossing of this area generally results in the excitation of the bodies’ eccentricities, leading to a collision or scattering instead of the capture into the co-orbital resonance.

Two processes that can form co-orbital exoplanets were proposed by Laughlin & Chambers (2002): either planet-planet gravitational scattering, or accretion in situ at the L₄/L₅ points of a primary. The assumptions made on the gas-disc density profile in the scattering scenario can either lead to systems with a high diversity of mass ratios (Cresswell & Nelson 2008, 2009), or to equal mass co-orbitals when a density jump is present (Giuppone et al. 2012). In their model, Cresswell & Nelson (2008) form co-orbitals in over 30% of the generated planetary systems, with very low inclinations and eccentricities (e <  0.02). In several hydrodynamical simulations of the formation of the outer part of the solar system by Crida (2009), Uranus and Neptune ended up in co-orbital configuration, while both were trapped in the same MMR with Saturn. In the in situ scenario, different studies yielded different upper limits to the mass that can form at the L₄/L₅ equilibrium point of a giant planet: Beaugé et al. (2007) obtained a maximum mass of ∼0.6 M_⊕, while Lyra et al. (2009) obtained 5 − 15 M_⊕ planets in the tadpole area of a Jupiter-mass primary.

The growth and evolution of co-orbitals have also been studied. For existing co-orbitals, Cresswell & Nelson (2009) found that gas accretion increases the mass difference between the planets, with the more massive of the two reaching Jovian mass while the starving one stays below 70 M_⊕. They also found that inward migration tends to slightly increase the amplitude of libration of the co-orbital, while remaining within the tadpole domain. The divergence from the resonance accelerates during late migrating stages with low gas friction, which may lead to instabilities. Another study from Pierens & Raymond (2014) shows that equal mass co-orbitals (from super-Earths to Saturns) are heavily disturbed during the gap-opening phase of their evolution. Rodríguez et al. (2013) have also shown that in some cases long-lasting tidal evolution may perturb equal mass close-in systems.

In this work we aim to better understand what causes the stability or instability of the co-orbital configuration during the protoplanetary disc phase. To do this we study the effects that planet migration, through interactions with a protoplanetary gas disc, has on the evolution of the co-orbital resonance angle. We also examine the evolution of the eccentricities and inclinations of the co-orbitals as the planets migrate throughout the disc. Both the type I (when a planet is fully embedded in the protoplanetary disc), and the type II (when the planet is massive enough to significantly perturb the disc, i.e. open a gap in the disc) migration regimes are studied.

This paper is laid out as follows. In Sect. 2, we describe the co-orbital dynamics in the absence of dissipation. In Sect. 3, we develop an integrable analytical model of the co-orbital resonance under a generic dissipation and mass change, in the coplanar-circular case. The application to type I migration is performed in Sect. 4, first using a simple analytical model for the disc torque, then comparing to the long-term evolution in an evolving protoplanetary disc. The result of population synthesis and the effect of resonant chains on the co-orbital configuration will be discussed in that section as well. We then estimate the forces that are actually applied on the coorbital by performing hydrodynamical simulation in Sect. 5. Finally, in Sect. 6, we discuss the stability of the co-orbital configuration in the direction of the eccentricities and inclinations. We then draw our conclusions in Sect. 7.

2. Dynamics of the co-orbital resonance in the non-dissipative case

In this section we describe the co-orbital motion of two planets of masses m₁ and m₂ around a central star of mass m₀ without any dissipative forces. For both planets we define: μ_j = 𝒢(m₀ + m_j), and β_j = m₀m_j/(m₀ + m_j), where 𝒢 is the gravitational constant. We note a_j the semi-major axis of planet j, e_j its eccentricity, and I_j its inclination. We use Poincaré astrocentric coordinates for both planets:

$$\begin{array}{cccc}\hfill {\mathrm{\Lambda }}_j& ={\beta }_j\sqrt{{\mu }_j{a}_j},\hfill & \hfill {\lambda }_j& ={\lambda }_j,\hfill \\ \hfill x_j& =\sqrt{{\mathrm{\Lambda }}_j}\sqrt{1-\sqrt{1-e_j^2}}{\mathrm{e}}^{i{\varpi }_j},\hfill & \hfill {\stackrel{\sim }{x}}_j& =-i{\overline{x}}_j,\hfill \\ \hfill y_j& =\sqrt{{\mathrm{\Lambda }}_j}\sqrt{\sqrt{1-e_j^2}(1-cosI)}{\mathrm{e}}^{i{\mathrm{\Omega }}_j},\hfill & \hfill {\stackrel{\sim }{y}}_j& =-i{\overline{y}}_j,\hfill \end{array}$$(1)

where λ_j, ϖ_j and Ω_j are its mean longitude, longitude of the pericenter, and ascending node of each planet, and ${\overline{x}}_j$ and ${\overline{y}}_j$ are the complex conjugates of x_j and y_j, respectively. The Hamiltonian of the system reads, in these coordinates:

$$\begin{array}{cc}\hfill H& =H_{\mathrm{K}}({\mathrm{\Lambda }}_1,{\mathrm{\Lambda }}_2)\hfill \\ \hfill & +H_{\mathrm{P}}({\lambda }_1,{\lambda }_2,{\mathrm{\Lambda }}_1,{\mathrm{\Lambda }}_2,x_1,x_2,{\stackrel{\sim }{x}}_1,{\stackrel{\sim }{x}}_2,y_1,y_2,{\stackrel{\sim }{y}}_1,{\stackrel{\sim }{y}}_2),\hfill \end{array}$$(2)

where H_K is the Keplerian component and H_P is the perturbative component due to planet-planet interactions taking into account both direct and indirect effects. The Keplerian component depends only on Λ_j:

$$\begin{array}{c}\hfill H_{\mathrm{K}}=-{\displaystyle \sum _{j=1}^2\left(\frac{{\mu }_j^2{\beta }_j^3}{2{\mathrm{\Lambda }}_j^2}\right),}\end{array}$$(3)

whereas the perturbative component depends on all twelve Poincaré coordinates. We do not need to express the explicit form of H_P at this point but it could be seen as an expansion of the x_j and y_j variables around 0. H_P can be obtained for example using the algorithm developed in Laskar & Robutel (1995).

As we study here the 1:1 mean motion resonance, we place ourselves in the neighbourhood of the exact Keplerian resonance defined by ${\dot{\lambda }}_1={\dot{\lambda }}_2$. The equations canonically associated with the Hamiltonian (2) hence read:

$$\begin{array}{c}\hfill \frac{\partial H_{\mathrm{K}}}{\partial {\mathrm{\Lambda }}_1}({\mathrm{\Lambda }}_1,{\mathrm{\Lambda }}_2)=\frac{\partial H_{\mathrm{K}}}{\partial {\mathrm{\Lambda }}_2}({\mathrm{\Lambda }}_1,{\mathrm{\Lambda }}_2).\end{array}$$(4)

We note ${\mathrm{\Lambda }}_1^0$ and ${\mathrm{\Lambda }}_2^0$ are the solutions of Eqs. (4) and (3). ${\mathrm{\Lambda }}_1^0$ and ${\mathrm{\Lambda }}_2^0$ are uniquely determined if we choose the exact Keplerian resonance which has the same total angular momentum than the studied orbit. At first order in e and I, the total angular momentum reads:

$$\begin{array}{c}\hfill L={\mathrm{\Lambda }}_1+{\mathrm{\Lambda }}_2={\mathrm{\Lambda }}_1^0+{\mathrm{\Lambda }}_2^0.\end{array}$$(5)

We can thus express ${\mathrm{\Lambda }}_1^0$ and ${\mathrm{\Lambda }}_2^0$ as a function of L:

$$\begin{array}{c}\hfill {\mathrm{\Lambda }}_{\mathit{jL}}^0=m_{j}L/(m_1+m_2)+\mathcal{O}(\epsilon ),\end{array}$$(6)

where ε = 𝒪((m₁ + m₂)/m₀). We can also define the average mean-motion η associated to the exact Keplerian resonance by, at order 0 in ε:

$$\begin{array}{c}\hfill {\eta }_{\mathrm{L}}=\frac{{\mu }_j^2{\beta }_j^3}{{({\mathrm{\Lambda }}_j^0)}^3}={\mu }_0^2{\left(\frac{m_1+m_2}{L}\right)}^3,\end{array}$$(7)

where μ₀ = 𝒢m₀. We also define the associated semi-major axis:

$$\begin{array}{c}\hfill {\overline{a}}_{\mathrm{\Gamma }}=\frac{L^2}{{\mu }_0{(m_1+m_2)}^2}·\end{array}$$(8)

2.1. Averaging the Hamiltonian near the co-orbital resonance

Since the mean motions n_j of the two bodies are close at any given time, the quantity ζ = λ₁ − λ₂ evolves slowly with respect to the longitudes. The Hamiltonian (2) hence possesses 3 time-scales: a fast one, associated with the mean motion η and the mean longitudes, a semi-fast one, associated with the resonant frequency $\nu =\mathcal{O}(\sqrt{\epsilon })$ and the libration of the resonant angle ζ, and a slow time-scale (called secular), which is associated with the orbital precession and the variables x_j, ${\stackrel{\sim }{x}}_j$, y_j and ${\stackrel{\sim }{y}}_j$. To emphasise the separation of these time-scales, we process to the following canonical change of variables (x_j and y_j remain unchanged):

$$\begin{array}{c}\hfill \left(\begin{array}{c}\zeta \\ {\zeta }_2\end{array}\right)=\left(\begin{array}{cc}1& -1\\ 0& 1\end{array}\right)\left(\begin{array}{c}{\lambda }_1\\ {\lambda }_2\end{array}\right),\left(\begin{array}{c}Z\\ Z_2\end{array}\right)=\left(\begin{array}{cc}1& 0\\ 1& 1\end{array}\right)\left(\begin{array}{c}{\mathrm{\Lambda }}_1\\ {\mathrm{\Lambda }}_2\end{array}\right).\end{array}$$(9)

The Hamiltonian now reads:

$$\begin{array}{cc}\hfill \mathcal{H}& ={\mathcal{H}}_{\mathrm{K}}(Z,Z_2)\hfill \\ \hfill & +{\mathcal{H}}_{\mathrm{P}}(\zeta ,{\zeta }_2,Z,Z_2,x_1,x_2,{\stackrel{\sim }{x}}_1,{\stackrel{\sim }{x}}_2,y_1,y_2,{\stackrel{\sim }{y}}_1,{\stackrel{\sim }{y}}_2).\hfill \end{array}$$(10)

The separation between the time-scales allows for the averaging over the rapid angle ζ₂. Following Robutel & Pousse (2013), Robutel et al. (2015), we obtain the Hamiltonian:

$$\begin{array}{c}\hfill \overline{\mathcal{H}}={\mathcal{H}}_{\mathrm{K}}(Z,Z_2)+{\overline{\mathcal{H}}}_{\mathrm{P}}(\zeta ,Z,Z_2,x_j,{\stackrel{\sim }{x}}_j,y_j,{\stackrel{\sim }{y}}_j).\end{array}$$(11)

2.2. Circular coplanar case

In the coplanar circular case, 1D models of the 1:1 mean-motion (co-orbital) can be obtained taking x_j = y_j = 0 and developing the Hamiltonian (11) at second order in $Z-{\mathrm{\Lambda }}_1^0$ and $Z_2-({\mathrm{\Lambda }}_1^0+{\mathrm{\Lambda }}_2^0)$ (Robutel et al. 2015). The equation canonically associated with that Hamiltonian can be rewritten as a 2nd order differential equation (Érdi 1977; Robutel et al. 2015):

$$\begin{array}{c}\hfill \ddot{\zeta }=-3{\eta }^2\frac{m_1+m_2}{m_0}(1-{(2-2cos\zeta )}^{-3/2})sin\zeta .\end{array}$$(12)

The phase space of Eq. (12) is shown in Fig. 1.

Out of the four fixed points of Eq. (12), the collision (ζ = 0°) is not in the validity domain of the equation and will be ignored. ζ = 180° is the hyperbolic (unstable) L₃ Lagrangian equilibrium, while ζ = ±60° are elliptic (stable) configurations, the L₄ and L₅ Lagrangian equilibria. Orbits that librate around these elliptic equilibria are called tadpole, or trojan (in reference to Jupiter’s trojan swarms). Examples of trojan orbits are shown in red in Fig. 1. The separatrix emanating from L₃ (black curve) delimits trojan orbits from horseshoe orbits (examples are shown in blue), for which the system undergoes large librations that encompasses L₃, L₄ and L₅.

As previously stated, the libration of the resonant angle ζ is slow with respect to the average mean-motion η. The fundamental libration frequency ν is proportional to $\sqrt{(m_1+m_2)/m_0}\eta$. In the neighbourhood of the L₄/L₅ equilibria (Charlier 1906):

$$\begin{array}{c}\hfill {\nu }_0=\sqrt{\frac{27}{4}}\sqrt{\frac{m_1+m_2}{m_0}}\eta .\end{array}$$(13)

2.3. Dynamics in the eccentric and inclined directions

In order to study the co-orbital dynamics for low eccentricities and inclinations, ℋ_P can be expanded in Taylor series in the neighbourhood of (x₁,x₂,y₁,y₂) = (0,0,0,0), at 2nd order in x_j, y_j. This expansion can be written as (Laskar 1989; Robutel & Pousse 2013):

$$\begin{array}{cc}\hfill \overline{\mathcal{H}}& ={\overline{\mathcal{H}}}_0(\zeta ,Z,Z_2)+{\overline{\mathcal{H}}}_x^{(2)}(\zeta ,Z,Z_2,x_j,{\stackrel{\sim }{x}}_j)\hfill \\ \hfill & +{\overline{\mathcal{H}}}_y^{(2)}(\zeta ,Z,Z_2,y_j,{\stackrel{\sim }{y}}_j),\hfill \end{array}$$(14)

where ${\overline{\mathcal{H}}}_0$ is given in Appendix A, and ${\mathcal{H}}_x^{(2)}$ and ${\mathcal{H}}_y^{(2)}$ are sums of quadratic monomials in (x_j,${\stackrel{\sim }{x}}_j$) and (y_j,${\stackrel{\sim }{y}}_j$), respectively.

From this form we learn two things: at low eccentricity and inclination, the dynamics of the variables x = (x₁, x₂) and y = (y₁, y₂) are decoupled, and the dynamics of the variables ζ and Z are independent of x_j and y_j at first order. In the conservative case, the variational equations of x and y are given by Robutel & Pousse (2013), Robutel et al. (2015):

$$\begin{array}{c}\hfill \dot{\mathit{x}}=M_x(\zeta )\mathit{x},\dot{\mathit{y}}=M_y(\zeta )\mathit{y},\end{array}$$(15)

with

$$\begin{array}{c}\hfill M_x(\zeta )=\left(\begin{array}{cc}\frac{A_x(\zeta )}{m_1}& \frac{{\overline{B}}_x(\zeta )}{\sqrt{m_1{m}_2}}\\ \frac{B_x(\zeta )}{\sqrt{m_1{m}_2}}& \frac{A_x(\zeta )}{m_2}\end{array}\right),M_y=\left(\begin{array}{cc}\frac{A_y(\zeta )}{m_1}& \frac{{\overline{B}}_y(\zeta )}{\sqrt{m_1{m}_2}}\\ \frac{B_y(\zeta )}{\sqrt{m_1{m}_2}}& \frac{A_y(\zeta )}{m_2}\end{array}\right),\end{array}$$(16)

where A_x, B_x, A_y, B_y are complex-valued functions of ζ and are given in Appendix A. $\overline{B}$ represents the complex conjugate of B.

2.3.1. The eccentric direction

The dynamics in the direction of the eccentricities is given by the system (15). Although the coefficients of the matrix M_x are functions of the resonant angle ζ, ζ does not depend on the eccentricities at first order. One can hence evaluate M_x at a fixed point of the (ζ,Z) variables (L₄ will be described here, results are equivalent for L₅), and study the dynamics of the x_j variable in the neighbourhood of the L₄ circular equilibria. The eigenvector of the matrix M_x(L₄) gives the direction of two remarkable families of quasi-periodic orbits that are represented in Fig. 2 (Giuppone et al. 2010; Robutel & Pousse 2013):

• The first eigenvector, paired with a null eigenvalue g₋ = 0, is tangent to the eccentric Lagrangian family (EL₄), for which e₁ = e₂ and Δϖ = ϖ₁ − ϖ₂ = ζ.

• The second eigenvector, paired with the eigenvalue g₊ = 27/8(m₁ + m₂)/m₀η, is tangent to the anti-Lagrangian family (AL₄). For low eccentricities, this family is tangent to m₁e₁ = m₂e₂ and Δϖ = ϖ₁ − ϖ₂ = ζ + π.

Description of the dynamics at larger eccentricities can be found for example in Nesvorný et al. (2002), Giuppone et al. (2010), and Leleu et al. (2018).

2.3.2. The inclined direction

The dynamics in the direction of inclination is given by the system (15). In opposition to the eccentric direction, we do not learn anything by evaluating M_y at the L₄ equilibria since all of its coefficients vanish for ζ = 60°. However, since the evolution of ζ is fast with respect to the secular evolution on the y_j variables, we can obtain an approximation of the secular dynamics in the direction of the inclination for a given trajectory by averaging the expression of this system over a period 2π/ν with respect to the time t.

We note that Im(B_y(ζ)) = −A_y(ζ) (see Appendix A). In addition, the real part of B_y(ζ) is proportional to the expression of $\ddot{\zeta }$ (Eq. (12)), which is the derivative of a periodic function of period 2π/ν. Its average value over 2π/ν is hence null. As a result, the averaged M_y can be written:

$$\begin{array}{c}\hfill {\overline{M}}_y=-i\frac{m_1{m}_2}{2m_0}\eta {\overline{A}}_y\left(\begin{array}{cc}\frac{1}{m_1}& \frac{-1}{\sqrt{m_1{m}_2}}\\ \frac{-1}{\sqrt{m_1{m}_2}}& \frac{1}{m_2}\end{array}\right),\end{array}$$(17)

where ${\overline{A}}_y$ is the averaged value of A_y(ζ) over a period 2π/ν, see Fig. 3. At ζ₀ = 60° (Lagrangian equilibrium), ${\overline{A}}_y=0$ and the system is degenerate. If ${\overline{A}}_y\ne 0$, we identify the two eigenvectors:

• The first eigenvector, paired with a null eigenvalue s₋ = 0, is tangent to the direction I₂ = I₁ and Ω₂ = Ω₁: the two planets orbit in the same plane, inclined by I₁ = I₂ with respect to the reference frame.

• The second eigenvector, paired with the eigenvalue $s_{+}=-i(m_1+m_2)\eta {\overline{A}}_y/(2m_0)$, is tangent to the direction m₁I₁ = m₂I₂ and Ω₂ = Ω₁ + π: the inclination of both co-orbitals is constant and their lines of nodes slowly precess at the frequency s₊.

We note that in the second direction the Oxy plane of the reference frame is perpendicular to the total angular momentum of the system (i.e. the Oxy plane is the invariant plane).

3. Stability of the Lagrangian equilibria L₄ and L₅ under external forces and mass change – the coplanar circular case

In this section we study the stability of the Lagrangian equilibria in the coplanar circular case under a generic dissipation, modelled by forces F₁ and F₂ applied on each planet. We assume that these forces are small with respect to the attraction by the central star. Using Gauss’ equations, the Poincaré variables are modified in the following way:

$$\begin{array}{cc}\hfill {\dot{\mathrm{\Lambda }}}_{j,d}& ={\mathrm{\Gamma }}_j=\frac{{\mathrm{\Lambda }}_j^2}{\mathcal{G}{m}_0{m}_j^2}{F}_{\mathrm{t},j},\hfill \\ \hfill {\dot{\lambda }}_{j,d}& =R_j=-\frac{2{\mathrm{\Lambda }}_j}{\mathcal{G}{m}_0{m}_j^2}{F}_{\mathrm{r},j}.\hfill \end{array}$$(18)

Where F_r, j is the radial force applied on the planet j, while Γ_j is the torque induced by the tangential force F_t, j. If the forces vary significantly over the orbital time-scale, they will also excite the eccentricities of the planets. For this work, we assume that variation of these forces over an orbital period is negligible, and that non-axisymmetric dissipative forces applied on the planets can be parametrised by Λ₁, Λ₂, and ζ = λ₁ − λ₂.

3.1. Model of the 1:1 MMR under dissipation and mass change

We start by developing an analytical model for the dynamics of the co-orbital configuration in the neighbourhood of the Lagrangian equilibria L₄ and L₅ in the dissipative case. To develop this model, we head back to the Hamiltonian transformations that we described in Sect. 2: we rewrite the equation of variation using the canonical variables Λ_j, λ_j. The equations of motion are given by the equation canonically associated with the Hamiltonian H, Eq. (2), to which we add the effect of the dissipation (Eq. (18)), and a slow, isotropic mass change for the planets parametrised by two constants ṁ₁ and ṁ₂.

We then perform a change of variables to uncouple the fast (i.e. associated to the mean motion) and semi-fast (i.e. resonant) degrees of freedom (Robutel et al. 2015):

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \left(\begin{array}{c}\zeta \\ \phi \end{array}\right)& =\left(\begin{array}{cc}1& -1\\ \frac{m_1}{m_1+m_2}& \frac{m_2}{m_1+m_2}\end{array}\right)\left(\begin{array}{c}{\lambda }_1\\ {\lambda }_2\end{array}\right),\hfill \\ \hfill \left(\begin{array}{c}\widehat{\mathrm{\Delta }}\\ L\end{array}\right)& =\left(\begin{array}{cc}\frac{m_2}{m_1+m_2}& -\frac{m_1}{m_1+m_2}\\ 1& 1\end{array}\right)\left(\begin{array}{c}{\mathrm{\Lambda }}_1\\ {\mathrm{\Lambda }}_2\end{array}\right).\hfill \end{array}\end{array}$$(19)

It should be noted that this change of co-ordinates is canonical in the conservative case, here we need to add the terms relative to the dissipative forces and mass changes. Averaging over the fast angle φ, we obtain the following system:

$$\begin{array}{cc}\hfill \dot{\widehat{\mathrm{\Delta }}}& ={\mu }_0^2\frac{m_1{m}_2{(m_1+m_2)}^2}{m_0{L}^2}(1-\frac{1}{\delta {(\zeta )}^3})sin\zeta \hfill \\ \hfill & +\frac{m_2{\mathrm{\Gamma }}_1-m_1{\mathrm{\Gamma }}_2}{(m_1+m_2)}+\frac{m_1^2{\dot{m}}_2+m_2^2{\dot{m}}_1}{m_1{m}_2(m_1+m_2)}\widehat{\mathrm{\Delta }},\hfill \\ \hfill \dot{\zeta }& =-3{\mu }_0^2\frac{{(m_1+m_2)}^5}{m_1{m}_2{L}^4}\widehat{\mathrm{\Delta }}+R_1-R_2,\hfill \\ \hfill \dot{L}& ={\mathrm{\Gamma }}_1+{\mathrm{\Gamma }}_2+\left(\frac{{\dot{m}}_1+{\dot{m}}_2}{m_1+m_2}\right)L,\hfill \\ \hfill \dot{\phi }& =f(\widehat{\mathrm{\Delta }},L,\zeta )+\frac{m_2{R}_2+m_1{R}_1}{m_1+m_2},\hfill \end{array}$$(20)

with f a polynomial function of $\widehat{\mathrm{\Delta }}$ and L, with trigonometric terms in ζ, and $\delta (\zeta )=\sqrt{2-2cos\zeta }$. Γ_j and R_j remain unchanged as we neglected their evolution over an orbital period. For the rest of this study, we focus on the evolution of the resonant degree of freedom ($\widehat{\mathrm{\Delta }}$,ζ).

3.2. Constant torques, radial forces and masses

At first we only consider the constant part of the torques and radial forces, Γ_j0 and R_j0, and constant masses (ṁ_j = 0). In this case, the evolution of the resonant degree of freedom ($\widehat{\mathrm{\Delta }}$,ζ) is given by the equations associated to the following Hamiltonian:

$$\begin{array}{cc}\hfill {\mathcal{H}}_{\mathrm{r}}& =\frac{m_1{m}_2{\eta }_{\mathrm{L}}L}{m_0(m_1+m_2)}(F(\zeta )-C_{\mathrm{\Gamma }}\zeta )\hfill \\ \hfill & -\frac{3}{2}{\eta }_{\mathrm{L}}\frac{{(m_1+m_2)}^2}{m_1{m}_{2}L}\widehat{\mathrm{\Delta }}(\widehat{\mathrm{\Delta }}-2{\widehat{\mathrm{\Delta }}}_{\mathrm{eq}}),\hfill \end{array}$$(21)

where F(ζ) = cosζ − (2 − 2cosζ)^−1/2,

$$\begin{array}{c}\hfill C_{\mathrm{\Gamma }}=(\frac{{\mathrm{\Gamma }}_{10}}{m_1}-\frac{{\mathrm{\Gamma }}_{20}}{m_2})\frac{m_0}{{\eta }_{\mathrm{L}}L},\end{array}$$(22)

and

$$\begin{array}{c}\hfill {\widehat{\mathrm{\Delta }}}_{\mathrm{eq}}=\frac{L{m}_1{m}_2(R_{10}-R_{20})}{3{\eta }_{\mathrm{L}}{(m_1+m_2)}^2}·\end{array}$$(23)

3.2.1. Asymmetry of the phase space

In the restricted case (m₁ ≪ m₀, m₂ = 0), the application of a constant torque on the co-orbitals result in a distortion of the phase space (Murray 1994; Sicardy & Dubois 2003). For a negative torque applied on the massive planet, it leads to a smaller tadpole domain for trailing massless particles than for leading ones, as the hyperbolic Lagrangian point L₃ gets closer to L₅, and further away from L₄. We study here the displacement of the perturbed circular coplanar Lagrangian equilibria for two massive bodies, as a function of the dimensionless quantity C_Γ, equivalent to the variable “α” in Sicardy & Dubois (2003). The position of these equilibria are obtained by solving the system:

$$\begin{array}{cc}\hfill \dot{\zeta }& =+\partial {\mathcal{H}}_{\mathrm{r}}/\partial \widehat{\mathrm{\Delta }}=0,\hfill \\ \hfill \dot{\widehat{\mathrm{\Delta }}}& =-\partial {\mathcal{H}}_{\mathrm{r}}/\partial \zeta =0.\hfill \end{array}$$(24)

For small enough dissipative forces (C_Γ ≪ 1), we can compute analytically the positions of the new equilibria in the neighbourhood of their position in the absence of dissipation, ζ = π/3, and $\widehat{\mathrm{\Delta }}=0$. In order to keep track of the relative size of the terms in Eq. (20), we introduce the small dimensionless tracer ε. We make the assumption that the perturbative terms C_Γ and $\sqrt{R_j/\eta }$ are of similar size with m_j/m₀, traced by ε. ε is just a tool to neglect second-order perturbative terms, and we can later take ε = 1 for numerical estimation of the variables. Using the implicit function theorem, we look for a shift of size ε in the value of ζ and $\widehat{\mathrm{\Delta }}$ with respect to the non-dissipative case. We hence replace $\widehat{\mathrm{\Delta }}$ by $\epsilon {\widehat{\mathrm{\Delta }}}_{L_4}$ and ζ by π/3 + εz in the system (24), and solve it. At lowest order in ε, the new L₄ equilibrium is located at:

$$\begin{array}{c}\hfill {\widehat{\mathrm{\Delta }}}_{L_4}={\epsilon }^3{\widehat{\mathrm{\Delta }}}_{\mathrm{eq}},z_{L_4}=-\frac{4}{9}\epsilon C_{\mathrm{\Gamma }}.\end{array}$$(25)

Similarly, we compute the evolution of the position of the L₃ equilibria by developing the system (24) in the neighbourhood of Δ = 0 and z = ζ − π = 0. Under the effect of the dissipative terms, the fixed point L₃ is shifted by:

$$\begin{array}{c}\hfill {\widehat{\mathrm{\Delta }}}_{L_3}={\epsilon }^3{\widehat{\mathrm{\Delta }}}_{\mathrm{eq}},z_{L_3}=\frac{8}{7}\epsilon C_{\mathrm{\Gamma }}.\end{array}$$(26)

We hence have:

$$\begin{array}{cc}\hfill {\zeta }_{L_3}-{\zeta }_{L_4}& =\frac{2\pi }{3}+\epsilon (z_{L_3}-z_{L_4})\hfill \\ \hfill & =\frac{2\pi }{3}+\frac{100}{63}\epsilon C_{\mathrm{\Gamma }}.\hfill \end{array}$$(27)

As a result, L₃ and L₄ get closer if $\frac{{\mathrm{\Gamma }}_{10}}{m_1}<\frac{{\mathrm{\Gamma }}_{20}}{m_2}$. If for example the torque per mass unit of the leading planet is lower that the torque per mass unit of the trailing planet, this results in a smaller trojan domain for the studied configuration.

For larger C_Γ, we solve the system (24) numerically. This system has three solutions for ζ ∈ [0 ° ,360 ° ] as long as |C_Γ|≲0.72. For larger absolute values of C_Γ, two of the three roots merge and vanish. The positions of the equilibria are shown in Fig. 4. For these 3 equilibria, the value of $\widehat{\mathrm{\Delta }}$ remain ${\widehat{\mathrm{\Delta }}}_{\mathrm{eq}}$. We also compute, for this value of $\widehat{\mathrm{\Delta }}$, the position of the separatrix emanating from L₃. To do so, we find the solutions of the equation ${\mathcal{H}}_{\mathrm{r}}({\widehat{\mathrm{\Delta }}}_{\mathrm{eq}},{\zeta }_{L_3})={\mathcal{H}}_{\mathrm{r}}({\widehat{\mathrm{\Delta }}}_{\mathrm{eq}},\zeta )$. The two separatrices were added as dashed lines to Fig. 4, and illustrate clearly the variation of the width of the trojan domain in the direction of ζ as a function of C_Γ, as the orbits librating around L₄ (resp. L₅) have to remain between the dashed and solid black lines. These results are in agreement with those of Sicardy & Dubois (2003), and generalise them to the case of two massive bodies.

3.2.2. Stability of the Lagrangian equilibria

We head back to small perturbations (C_Γ ≪ 1, $\sqrt{R_j/\eta }\ll 1$). The torques applied on each planet slowly change the total angular momentum of the system. Here we describe how the evolution of L changes the orbit of the co-orbitals on a long time-scale with respect to the resonant motion.

Let us consider a pair of co-orbitals on a trajectory librating in the neighbourhood of the L₄ equilibria. We use z = ζ − ζ_L4 and $d\widehat{\mathrm{\Delta }}=\widehat{\mathrm{\Delta }}-{\widehat{\mathrm{\Delta }}}_{\mathrm{eq}}$ to describe the trajectory. On the resonant time-scale, this trajectory follow a level curve of the Hamiltonian (21). This level curve ℒ can be parametrised by z_ℒ >  0, with the energy of the curve being H_r(Δ_eq, ζ_L4 + z_ℒ). For z_ℒ = 0, ℒ(z_ℒ) goes through $d\widehat{\mathrm{\Delta }}={\widehat{\mathrm{\Delta }}}_{\mathcal{L}}$, estimated by solving ${\mathcal{H}}_{\mathrm{r}}({\widehat{\mathrm{\Delta }}}_{\mathrm{eq}},{\zeta }_{L_4}+z_{\mathcal{L}})={\mathcal{H}}_{\mathrm{r}}({\widehat{\mathrm{\Delta }}}_{\mathrm{eq}}+{\widehat{\mathrm{\Delta }}}_{\mathcal{L}},{\zeta }_{L_4})$. We obtain:

$$\begin{array}{c}\hfill {\widehat{\mathrm{\Delta }}}_{1\mathcal{L}}=\frac{\sqrt{3}}{2}{\epsilon }^{3/2}W{z}_{\mathcal{L}}+\mathcal{O}({\epsilon }^{5/2}{z}_{\mathcal{L}}),\end{array}$$(28)

where

$$\begin{array}{c}\hfill W=\frac{m_1{m}_{2}L}{\sqrt{m_0}{(m_1+m_2)}^{3/2}}·\end{array}$$(29)

As L slowly evolves, so does W and H_r, and hence the level curve followed by the trajectory. However, the area enclosed by this level curve is an adiabatic invariant (Henrard 1982). As a result, a slow decrease of W leads the trajectory to follow level curves of larger and larger z_ℒ, while ${\widehat{\mathrm{\Delta }}}_{1\mathcal{L}}$ decreases, see the left panel of Fig. 5 for an example. We hence need to give a more precise definition of the stability we want to consider. As the position of the fixed points and separatrix are scale-free in the ζ direction, see Fig. 4, we consider a trajectory to be converging toward the Lagrangian equilibria if z_ℒ decreases, as the trajectory is getting further away from the boundary of the tadpole orbits.

We hence normalise the variable $\widehat{\mathrm{\Delta }}$ by W, and call this new variable Δ:

$$\begin{array}{c}\hfill \mathrm{\Delta }=\frac{\widehat{\mathrm{\Delta }}}{W}·\end{array}$$(30)

By doing so, we lose the Hamiltonian formulation of the problem as we introduce a dissipative term in the equations of variation, but we explicit the effect of the change of total angular momentum on the newly-defined stability of the system: despite the small dissipative term, on the resonant time-scale a trajectory will remain close to a level curve of the Hamiltonian part of the system. For these new level curves, Δ_ℒ only depends on z_ℒ. In these new variables, the stability is hence defined as a convergence toward L₄ in both the ζ and Δ directions, while divergence from L₄ happens in both ζ and Δ directions as well, as we can see in the right panel of Fig. 5. The evolution of Δ is given by:

$$\begin{array}{c}\hfill \dot{\mathrm{\Delta }}=\frac{1}{W}\dot{\widehat{\mathrm{\Delta }}}-\frac{{\dot{W}}_{\mathrm{c}}}{W}\mathrm{\Delta },\end{array}$$(31)

where

$$\begin{array}{c}\hfill \frac{{\dot{W}}_{\mathrm{c}}}{W}=\frac{{\mathrm{\Gamma }}_{10}+{\mathrm{\Gamma }}_{20}}{L}·\end{array}$$(32)

where Ẇ_c is the term of Ẇ coming from the constant torques applied on the planets. The resonant part of the system of equations of variation (20) becomes:

$$\begin{array}{cc}\hfill \dot{\mathrm{\Delta }}& ={\eta }_{\mathrm{L}}\sqrt{\frac{m_1+m_2}{m_0}}(1-\frac{1}{\delta {(\zeta )}^3})sin\zeta \hfill \\ \hfill & -\frac{{\dot{W}}_{\mathrm{c}}}{W}\mathrm{\Delta }+\frac{m_2{\mathrm{\Gamma }}_{10}-m_1{\mathrm{\Gamma }}_{20}}{(m_1+m_2)W},\hfill \\ \hfill \dot{\zeta }& =-3{\eta }_{\mathrm{L}}\sqrt{\frac{m_1+m_2}{m_0}}\mathrm{\Delta }+R_{10}-R_{20}.\hfill \end{array}$$(33)

As this is our main set of variables, we give the expression of the variable Δ with respect to the orbital elements:

$$\begin{array}{cc}\hfill \mathrm{\Delta }& \equiv \frac{\sqrt{m_0}\sqrt{m_1+m_2}(m_2{\mathrm{\Lambda }}_1-m_1{\mathrm{\Lambda }}_2)}{m_1{m}_{2}L}\hfill \\ \hfill & \equiv \sqrt{{\mu }_0{m}_0}\frac{\sqrt{m_1+m_2}}{L}(\sqrt{a_1}-\sqrt{a_2}).\hfill \end{array}$$(34)

while reciprocally, the circular angular momentum of each planet reads:

$$\begin{array}{cc}\hfill {\mathrm{\Lambda }}_1& =\frac{m_1}{m_1+m_2}L(1+\frac{m_2}{\sqrt{m_0(m_1+m_2)}}\mathrm{\Delta }),\hfill \\ \hfill {\mathrm{\Lambda }}_2& =\frac{m_2}{m_1+m_2}L(1-\frac{m_1}{\sqrt{m_0(m_1+m_2)}}\mathrm{\Delta }).\hfill \end{array}$$(35)

We now (and for the rest of this paper) study the stability of the new L₄ equilibria (ζ ∈ [0, 180 ° ]), while the stability of L₅ can be studied by swapping the indices of the planets. The L₄ equilibria is located in:

$$\begin{array}{c}\hfill \{\begin{array}{c}{\mathrm{\Delta }}_{L_4}=\frac{1}{3}{\epsilon }^{3/2}\frac{1}{{\eta }_{\mathrm{L}}}\sqrt{\frac{m_0}{m_1+m_2}}(R_{10}-R_{20}),\hfill \\ {\zeta }_{L_4}=\frac{\pi }{3}-\frac{4}{9}\epsilon C_{\mathrm{\Gamma }}.\hfill \end{array}\end{array}$$(36)

To do so, we linearise the resonant part of the system (20) in the neighbourhood of Δ_L4 and ζ_L4:

$$\begin{array}{c}\hfill \left(\begin{array}{c}\dot{\mathrm{\Delta }}\prime \\ z\prime \end{array}\right)=J_{L_4}\left(\begin{array}{c}\mathrm{\Delta }\prime \\ z\prime \end{array}\right),\end{array}$$(37)

where Δ′=Δ − Δ_L4, z′ = z − z_L4, and J₄ is the Jacobian matrix of the system (20) computed at the equilibrium (36). We make a final change of coordinates that diagonalises the system (37). In this new set of variables (z₁,z₂), the equations of variation reads:

$$\begin{array}{c}\hfill \left(\begin{array}{c}{\dot{z}}_1\\ z_2\end{array}\right)=\left(\begin{array}{cc}{u}_{\mathrm{c}}-i\nu & 0\\ 0& u_{\mathrm{c}}+i\nu \end{array}\right)\left(\begin{array}{c}{z}_1\\ z_2\end{array}\right)\end{array}$$(38)

where:

$$\begin{array}{c}\hfill \nu ={\epsilon }^{1/2}{\nu }_0+\mathcal{O}({\epsilon }^{3/2}),\end{array}$$(39)

and

$$\begin{array}{c}\hfill 2u_{\mathrm{c}}=-{\epsilon }^2\frac{{\dot{W}}_{\mathrm{c}}}{W}·\end{array}$$(40)

We hence obtain a modification of the classical resonant frequency ν₀ in the neighbourhood of L₄ (39), plus a hyperbolic term given in (40). As u_c and ν are not constant due to the evolution of the masses and the total angular momentum L, we will study the stability of “partial” equilibria (see for example Vorotnikov 2002) by dividing the variables into two groups: the variables with respect to which the stability is investigated (Δ,ζ), and the remaining variables L, m₁ and m₂. The linearised system with a diagonal matrix Eq. (38) allows for a trivial application of results on the stability of the partial equilibrium z₁ = z₂ = 0 (see Appendix B): The partial equilibrium (36) is stable if u_c is negative or null. From now on, we also make the assumption that a positive value induces a divergence from this equilibrium.

As previously stated, Ẇ_c/W represents the change of the size of the resonance as L evolves. If Ẇ_c/W < 0, the width of the resonance in the previous variable $\widehat{\mathrm{\Delta }}$ is decreasing, leading to a slow increase of z_l to retain a constant area enclosed by the level curve of the trajectory. As the system increases its amplitude of libration of the resonant angle, and gets closer to the separatrix of the tadpole domain, we consider the system to be diverging from the equilibrium. Looking at the expression of Ẇ_c/W, Eq. (32), it appears explicitly that inward migration (negative torques) have a destabilising effect on the co-orbital resonance, while outward migration tends to stabilise it.

3.3. Stability criteria for the Lagrangian configuration L₄ under non-constant forces and mass change

In this section we study the stability of the Lagrangian equilibria L₄ assuming that the variation of the dissipative forces can be parametrised by ζ and Δ, and that other variations are slow enough to be considered as constant with respect to the resonant time-scale. We also consider a slow, isotropic mass change for both masses. As in Sect. 3.2, we make the assumption that the perturbative terms C_Γ, $\sqrt{{\dot{m}}_j/m_j}$ and $\sqrt{R_j/\eta }$ are of similar size with m_j/m₀, traced by ε. As we did in the previous section, we normalise $\widehat{\mathrm{\Delta }}$ by W. As we now consider ṁ_j ≠ 0, we have:

$$\begin{array}{c}\hfill \frac{\dot{W}}{W}=\frac{{\mathrm{\Gamma }}_{01}+{\mathrm{\Gamma }}_{02}}{L}+\frac{\frac{{\dot{m}}_1}{m_1}(m_1+2m_2)+\frac{{\dot{m}}_2}{m_2}(m_2+2m_1)}{2(m_1+m_2)}·\end{array}$$(41)

The equations of variations (20) becomes:

$$\begin{array}{cc}\hfill \dot{\mathrm{\Delta }}& ={\eta }_{\mathrm{L}}\sqrt{\frac{m_1+m_2}{m_0}}(1-\frac{1}{\delta {(\zeta )}^3})sin\zeta \hfill \\ \hfill & -\frac{\dot{W}}{W}\mathrm{\Delta }+\frac{m_2{\mathrm{\Gamma }}_1-m_1{\mathrm{\Gamma }}_2}{(m_1+m_2)W}+\frac{m_1^2{\dot{m}}_2+m_2^2{\dot{m}}_1}{m_1{m}_2(m_1+m_2)}\mathrm{\Delta },\hfill \\ \hfill \dot{\zeta }& =-3{\eta }_{\mathrm{L}}\sqrt{\frac{m_1+m_2}{m_0}}\mathrm{\Delta }+R_1-R_2.\hfill \end{array}$$(42)

For the variables $\mathrm{\Delta }=\widehat{\mathrm{\Delta }}/W$, ζ. Then we linearise these perturbations in the neighbourhood of (Δ_L4,ζ_L4) by using the reduced variables Δ′=Δ − Δ_L4, z = ζ − ζ_L4. Γ_j and R_j become:

$$\begin{array}{cc}\hfill {\mathrm{\Gamma }}_j& ={\mathrm{\Gamma }}_{j0}+{\mathrm{\Gamma }}_{j\mathrm{\Delta }}\mathrm{\Delta }\prime +{\mathrm{\Gamma }}_{j\zeta }z,\hfill \\ \hfill R_j& =R_{j0}+R_{j\mathrm{\Delta }}\mathrm{\Delta }\prime +R_{j\zeta }z,\hfill \end{array}$$(43)

where Γ_jΔ = ∂Γ_j/∂Δ′, Γ_jζ = ∂Γ_j/∂z, and similar expressions for R_j. Despite the addition of these new terms in the equations of variation, the dominant terms of Δ_L4 and ζ_L4 remains those computed in the case of constant forces, given in (36).

We hence linearise the rest of the system (42) in the neighbourhood of (Δ_L4,ζ_L4), modelling the dissipative forces by the expression (43). Then, as in Sect. 3.2.2, we diagonalise the linearised system. Since the masses are not constant in this section, this transformation adds additional terms proportional to ṁ₁ and ṁ₂. However, as these terms are of size 𝒪(ε³), we can neglect them when we compute the eigenvalues of the linearised system. The dominant terms of the imaginary part remain:

$$\begin{array}{c}\hfill \nu ={\epsilon }^{1/2}{\nu }_0+\mathcal{O}({\epsilon }^{3/2}).\end{array}$$(44)

However, the real part of the eigenvalue gets new terms:

$$\begin{array}{cc}\hfill 2u& ={\epsilon }^2[R_{1\zeta }-R_{2\zeta }+\frac{m_2{\mathrm{\Gamma }}_{1\mathrm{\Delta }}-m_1{\mathrm{\Gamma }}_{2\mathrm{\Delta }}}{(m_1+m_2)W}\hfill \\ \hfill & -\frac{{\mathrm{\Gamma }}_{01}+{\mathrm{\Gamma }}_{02}}{W}-\frac{1}{2}\frac{{\dot{m}}_1+{\dot{m}}_2}{m_1+m_2}].\hfill \end{array}$$(45)

We hence obtain a small modification of the classical resonant frequency ν₀ in the neighbourhood of L₄ (44), plus a hyperbolic term given in Eq. (45). We note that at lowest order in ε the equilibrium point remains elliptic because of our assumptions on the relative size of the dissipative terms with respect to m_j/m₀. As in Sect. 3.2.2, we can deduce the stability of the system by estimating the sign of u: a negative value of u induces a convergence of the system toward the exact resonance, while a positive value led to a divergence.

The term-by-term physical interpretation of this criterion (the sign of u, Eq. (45)) is straightforward: if for example R_1ζ is the only non-zero term of u and using Eqs. (18) and (43), we have:

$$\begin{array}{c}\hfill \dot{z}=R_{1\zeta }z.\end{array}$$(46)

If R_1ζ <  0, z converges toward 0, hence ζ converges toward ζ_L4. Similarly, a negative Γ_1Δ implies that the planet “1” migrates inward faster when Δ >  0 (i.e. when a₁ >  a₂), which also pushes the planet towards the exact resonance.

The terms in Γ_j0 and ṁ_j were already described in the previous section, and take into account the evolution of the width of the resonance as the masses and total angular momentum of the configuration slowly evolve.

In the following sections, we apply the stability criterion u in the case of planets evolving in a protoplanetary disc.

4. Axi-symmetric dissipative forces: 1D protoplanetary discs

Gravitational interactions between the planets and their parent disc impact on the planets’ orbital parameters, typically causing them to migrate, either inwards towards the star or outwards away from the star (Goldreich & Tremaine 1979; Artymowicz 1993; Papaloizou & Larwood 2000). Planet eccentricities and inclinations are also affected by the interactions with the disc, typically causing them to be damped, forcing the planets to orbit their parent star on coplanar circular orbits (Cresswell & Nelson 2006; Bitsch & Kley 2010). The orbital evolution of the planets is the result of various torque components from the disc acting on them, such as the Lindblad torque, corotation torque, and horseshoe drag (Baruteau et al. 2014). For a single planet in a protoplanetary disc, two regimes are usually considered: As long as a planet is not massive enough to perturb the disc considerably, the planet is called to be in type I migration. As it is growing, it starts to open a gap around its orbit, and once this gap is deep enough, it migrates in type II regime. In this section, we assume that the perturbation of the disc by each planet is negligible, and hence that usual type I migration formulae can be applied on each planet individually.

We also assume that the unperturbed disc is axi-symmetric, and we study the stability of the coplanar circular co-orbital resonance in this case. Under this assumption, only the tangential forces (or torque) will play a role in the stability of the configuration (see the expression of the criterion u, Eq. (45)). In the literature, these torques are often modelled by migration time-scales τ_aj, used as prescriptions for the evolution of the semi-major axes: ${\dot{a}}_{j,d}=-a_j/{\tau }_{\mathit{aj}}$, which implies a dissipative term for the evolution of the angular momentum of the planet of the form: ${\dot{\mathrm{\Lambda }}}_{j,d}=-{\mathrm{\Lambda }}_j/(2{\tau }_{\mathit{aj}})$.

4.1. Analytic model of type 1 migration

Tanaka et al. (2002) and Tanaka & Ward (2004) derived a linear model of the wave excitation in three-dimensional isothermal discs to obtain an analytical model for the torques induced by the Lindblad and corotation resonances. Following their notation, we consider a gaseous disc parametrised by its aspect ratio h and its surface density Σ(r), such that:

$$\begin{array}{c}\hfill h=H/r=h_0{r}^f\end{array}$$(47)

where H is the scale height, h₀ is the aspect ratio at 1 au and f is the flaring index, and

$$\begin{array}{c}\hfill \mathrm{\Sigma }(r)={\mathrm{\Sigma }}_0{r}^{-\alpha }\end{array}$$(48)

where Σ₀ is the surface density at 1 au and α parametrises the slope of the surface density. At first order in the inclination and eccentricity, we have:

$$\begin{array}{c}\hfill {\tau }_{\mathit{aj}}=\frac{{\tau }_{\mathit{wj}}}{2.7+1.1\alpha }{h}^{-2}\end{array}$$(49)

for the evolution of the semi-major axis (Tanaka et al. 2002), with

$$\begin{array}{c}\hfill {\tau }_{\mathit{wj}}={\left(\frac{m_j}{m_0}\right)}^{-1}\frac{m_0}{\mathrm{\Sigma }(r_j)r_j^2}{h}^4{\mathrm{\Omega }}_j^{-1},\end{array}$$(50)

where ${\mathrm{\Omega }}_p={(\mathcal{G}{m}_0/r_j^3)}^{1/2}$ is the Keplerian angular velocity of the planet. Considering Eqs. (47) and (48), and neglecting the effect of the eccentricity on the star-planet distance, τ reads:

$$\begin{array}{c}\hfill {\tau }_{\mathit{wj}}=\frac{m_0}{m_j}\sqrt{\frac{m_0}{\mathcal{G}}}\frac{h_0^4}{{\mathrm{\Sigma }}_0}{a}_j^{\alpha +4f-1/2}.\end{array}$$(51)

At this stage, it becomes clear that taking τ_aj as a constant is a strong assumption that would force the semi-major axis to behave as an exponential decay. We hence introduce the variable K_j that parametrises the local slope of the migration:

$$\begin{array}{c}\hfill {\tau }_{\mathit{aj}}={\tau }_j{a}_j^{K_j}.\end{array}$$(52)

The effect of the disc on the angular momentum of each planet is hence given by:

$$\begin{array}{c}\hfill {\dot{\mathrm{\Lambda }}}_{j,d}=-m_j^{2K_j}{\mu }_0^{K_j}\frac{{\mathrm{\Lambda }}_j^{1-2K_j}}{2{\tau }_j}·\end{array}$$(53)

Expanding the torque applied on each planet at first order in Δ, ${\dot{\mathrm{\Lambda }}}_{j,d}={\mathrm{\Gamma }}_{j0}+{\mathrm{\Gamma }}_{j\mathrm{\Delta }}\mathrm{\Delta }$, we have:

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{j0}=-\frac{m_j}{m_1+m_2}\frac{L}{2{\overline{a}}_{\mathrm{L}}^{K_j}{\tau }_j},\end{array}$$(54)

and:

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{j\mathrm{\Delta }}={(-1)}^j\frac{(1-2K_j)L}{2{\overline{a}}_{\mathrm{L}}^{K_j}{\tau }_j}·\end{array}$$(55)

The effect of dissipative forces modelled by τ_aj (Eq. (52)) on the co-orbital configuration can be deduced from the expressions (54) and (55) and the criterion u, Eq. (45). The Lagrangian point is attractive if:

$$\begin{array}{c}\hfill ((2K_2-1)m_1+m_2){\tau }_1+((2K_1-1)m_2+m_1){\tau }_2<0.\end{array}$$(56)

Figure 6 represents the K values for which the equilibrium becomes repulsive (above the line of a given mass ratio) as a function of τ₁/τ₂ ≈ τ_a1/τ_a2, in the special case K = K₁ = K₂. These results apply to the neighbourhood of the L₄ or L₅ equilibria, hence to tadpole orbits with small amplitudes of libration, as they were derived from the linearisation of the system (42) in the neighbourhood of L₄. We then check the validity of this criterion for different amplitudes of libration in the case m₁ = 10m₂. We integrate the equations of the 3-body problem using the variable-step integrator DOPRI (Dormand & Prince 1980). In addition, the migration of the semi-major axis is modelled using

$$\begin{array}{c}\hfill {\ddot{\mathit{r}}}_j=-{\dot{\mathit{r}}}_j{a}_j^{-K}/(2{\tau }_j),\end{array}$$(57)

where r_j is the position of the planet j with respect to the star. The tests were made using m₀ = 1, m₁ = 5 × 10⁻⁵, m₂ = 5 × 10⁻⁶, taking as initial conditions e_j = I_j = 0, a₁ = a₂ = 1 au and ζ₀ = 58°, 50°, 40°, 30° and 20°. For each initial value of ζ, a grid of cases were integrated for different values of K and τ_w1/τ_w2, using τ_w1 = 1/m₁. The grey squares in Fig. 6 represent the threshold values of K for which the configurations change from converging (below the curve) to diverging (above it) for ζ₀ = 20° (horseshoe configuration). All other initial amplitudes of libration (initial values of ζ₀) gave very similar results. It implies that, at least for the chosen masses, the attraction of the exact resonance does not depend significantly on the amplitude of libration in the case of axi-symmetric dissipative forces.

In the case of the torque from Tanaka et al. (2002), we have

$$\begin{array}{c}\hfill {\tau }_{\mathit{aj}}={\tau }_j{a}_j^K=\frac{1}{2.7+1.1\alpha }\frac{m_0}{m_j}\sqrt{\frac{m_0}{\mathcal{G}}}\frac{h_0^2}{{\mathrm{\Sigma }}_0}{a}_j^{\alpha +2f-1/2}.\end{array}$$(58)

Which hence verifies the relation K₁ = K₂ = α + 2f − 1/2, in addition to τ_a1/τ_a2 = m₂/m₁. Applying these additional constraints, the stability of the Lagrangian points depends only on m₂/m₁ and the parameter K. The critical value for K is thus:

$$\begin{array}{c}\hfill K_0=-\frac{{(\frac{m_2}{m_1}-1)}^2}{4\frac{m_2}{m_1}}\end{array}$$(59)

which is represented in Fig. 7.

4.2. Stability in an evolving protoplanetary disc

The section above describes the evolution of the co-orbital resonance in a static environment, where the migration time-scales and disc parameters are assumed to be constant. However, protoplanetary discs do not remain static, so therefore it is important to determine how the co-orbital resonance evolves in a more global, ever-changing environment. Thus we now examine the behaviour of the resonance as the protoplanetary disc evolves on Myr time-scales, and also as a planet migrates from one region of a disc to another where the value of K can differ significantly. We use the disc model presented in Coleman & Nelson (2016b), where the standard diffusion equation for a 1D viscous α-disc model is solved (Shakura & Sunyaev 1973). Temperatures are calculated by balancing viscous heating and stellar irradiation against blackbody cooling. We use the torque formulae from Paardekooper et al. (2010, 2011) to calculate type I migration rates due to Lindblad and corotation torques acting on a planet. The Lindblad torque emerges when an embedded planet perturbs the local disc material, forming spiral density waves that are launched at the Lindblad resonances in the disc. Corotation torques arise from both local entropy and vortensity gradients in the disc, and their possible saturation is included in these simulations. The influence of eccentricity and inclination on the migration torques and the damping of eccentricities and inclinations are also included (Cresswell & Nelson 2008; Fendyke & Nelson 2014). These torques exchange angular momentum between the planet and the gas disc, and depending on their strength and direction can result in a torque being exerted on the planet, either inwards or outwards.

Figure 8 shows the gas surface density (left panel), disc aspect ratio h (middle panel) and the calculated value for K (right panel) at 4 different times throughout the disc lifetime. The blue line shows the profiles at the beginning of the disc lifetime, with the red, yellow and purple lines showing the profiles at 1, 2 and 3 Myr respectively. The disc lifetime here was ∼3.6 Myr. K = α + 2f − 1/2, as defined in Sect. 4.1.

We then examine the evolution of a pair of co-orbital planets that are undergoing type I migration in the protoplanetary disc represented in Fig. 8. The masses of the planets are 10 and 5 M_⊕ for the primary and secondary respectively, giving a mass ratio of 2. The top panel of Fig. 9 shows the evolution of semimajor axis as the planets migrate through the disc, showing that they remain in co-orbital resonance, even when migrating from 10 au down to the inner edge of the disc. The bottom panel shows the values for the surface density gradient α (blue line), the aspect ratio gradient f (red line), and the corresponding value of K (yellow line) at the planet’s location. As the planets originate in the outer irradiation dominated region of the disc, K has a value ∼1.75. But as the planets migrate in closer to the central star, they enter the viscous dominated region of the disc where the transition in opacities significantly impacts α and f, causing K to drop to below 0 for a time before rising back to just above 0. As the planets near the inner edge of the disc, K settles to just below 0. Whilst the planets are migrating through the different regions of the disc, they will either be converging towards the resonance or diverging away from the resonance, depending on the local disc conditions. Figure 10 shows the amplitude of libration (top panel) and the corresponding K′=K − K₀ for the two migrating planets in Fig. 9.4pt

When looking at Fig. 10, it can be seen that the amplitude of libration is increasing at the start of the simulation whilst K′∼2. As the planets migrate into the inner regions of the disc, K′ drops to begin fluctuating around 0. The first dashed line shows where the amplitude of libration begins to converge, and this lines up just after K′ reaches negative values. This is expected from the analytical model in Sect. 4.1. However, the analytical model uses the simplified type I migration torque formulae from Tanaka et al. (2002), whereas the evolving disc model here uses torque formulae from Paardekooper et al. (2010, 2011) that includes a more accurate treatment of the corotation torque. Despite these differences, it is interesting to see that the behaviour of the co-orbital resonance roughly matches what is expected from the Tanaka formulae, i.e. converging when K′ is negative and diverging when K′ is positive.

The planets migrating in Figs. 9 and 10 migrated until they reached the disc inner edge, close to the central star. Whilst doing so, they moved away from the co-orbital resonance, eventually breaking out of the resonance. However global simulations of planet formation have shown that planets in co-orbital configurations cover a wide range of semimajor axes at the end of the disc lifetime (Coleman & Nelson 2016a,b).

We now consider the case where the mass of a planet in the co-orbital resonance is changing over time in addition to its migration in the disc, meaning that it is accreting gas or planetesimals. As show in Sect. 3.3, both slow mass accretion and migration in the disc changes the width of the resonance in the way that leads to a divergence from the equilibrium (if the size of the resonance decreases due to inward migration), or to a convergence towards it (mass accretion or outward migration).

Figure 11 shows the relative change in mass against the relative change in angular momentum for a number of evolving planets at different locations of the protoplanetary disc and of different masses. We note that these planets are not in co-orbital configuration but are used to probe the torque felt and their accretion rate during their evolution in the disc. The planets shown in Fig. 11 were placed in a nominal protoplanetary disc similar to those shown in Coleman & Nelson (2016b). The black line shows where the sum relative changes in mass and angular momentum equate to zero, indicating no change in the stability of the co-orbital resonance. Planets that are evolving above this line will be accreting mass at a faster rate than they are migrating inwards and as such will be introducing a stabilising effect to their co-orbital region. For planets evolving below the line, then the opposite will occur, and they will diverge from the equilibrium. It is interesting to note, but not shown in the figure, that if a planet is undergoing outward migration, then the co-orbital resonance will be always be stabilising, so long as the planet is not losing mass at a significant rate. Looking at the regimes that specific planets operate in, we see that when a planet is small and growing through planetesimal or pebble accretion (yellow line), it is typically above the black line, since the planet is increasing in mass faster than it is migrating, stabilising the co-orbital region. For more massive planets, of mass between 10–20 M_⊕, that are accreting few pebbles and/or planetesimals, but are accreting gas slowly (red line), they will tend to sit below the black line, destabilising the resonance. This is due to them migrating faster than they are accreting, however this is also highly dependant on the local disc profiles, which could allow planets to become trapped and migrate inwards slowly, reducing the magnitude of the destabilisation. For example, the red line shows the evolution of a ∼15 M_⊕ planet initially orbiting at 5 au. The location of the planet in the protoplanetary disc, will also affect how close they appear to the black line in Fig. 11, since the migration torques are dependant on the local disc conditions, for example the 15 M_⊕ planet shown by the red line would be shifted to the right of the plot if it was initially orbiting at 20 au, instead of 5 au as is shown in the plot. For planets that grow into giant planets through runaway gas accretion (blue and purple lines), initially they are in a regime of significant stabilisation as they accrete gas extremely quickly. Once the runaway gas accretion phase ends, they accrete gas at a slower rate, but migrate at a similar rate, reducing the effects of the stabilisation, before ultimately making it destabilising. Again, given the location in the disc that these planets occupy, this will mainly affect the rate of change of angular momentum, possibly making the co-orbital region stabilise or destabilise at a faster rate. For the giant planet that survives migration (purple line), meaning it does not migrate into the central star, the period of slow type-II migration (at the bottom of Fig. 11) is destabilising for the co-orbital region. This is again due to the planet migrating at a faster relative rate than it is able to accrete gas.

Both mass accretion and migration hence have significant impacts on the evolution of the co-orbital resonance in the disc and have to be taken into account. In Sect. 5, we estimate how the perturbations of the disc induced by the pair of planets perturbs the torque that is applied to each of them.

4.3. Population outcome

While Sect. 4.2 examined the evolution of the co-orbital resonance in an evolving protoplanetary disc, it is interesting to see how often these co-orbital resonances actually occur in a much larger suite of simulations. To do this, we searched for co-orbital resonances in a recent population of simulations which studied planet formation around low mass stars. The simulations initially began with a realistic range of initial conditions, such as disc mass and solid mass, and used a similar disc model to that described in Sect. 4.2, but which had been adapted to being suitable for a low mass star (see Sect. 3 of Coleman et al. 2017). The main aim of these simulations was to examine the formation of planetary systems, through pebble or planetesimal accretion, around low-mass stars of mass 0.1 M_⊙, similar to Trappist-1 and Proxima Centauri (Coleman et al. 2019). Initially in these simulations, a number of low mass planetary embryos (m_p <  0.1 M_⊕) would be scattered throughout the disc, and would be able to either accrete pebbles or planetesimals, and also undergo type I migration. As the planets migrate they become trapped in resonant chains, typically involving first-order resonances, but could also become trapped in co-orbital resonances. As the systems evolve, the planets migrate to their final locations, sometimes maintaining their resonant chains and also their co-orbital configurations. The planetary systems are integrated for 3–5 million years after disc dispersal to allow the planetary systems to continue to evolve in an undamped environment.

We analysed the outcome of the 880 systems generated in the pebble accretion scenario and found co-orbitals in ≈12% of the final systems. Figure 12 displays all of these co-orbitals as a function of their semi-major axis and amplitude of libration at the end of the simulation. As expected, co-orbitals that migrated on their own and survived until the end of the disc lifetime tend to have a large amplitude of libration. However, the subgroup of co-orbitals that were trapped into a resonant chain with other planets seemed to be able to migrate close to the inner edge of the disc while retaining a small amplitude of libration.

Indeed, we show in Appendix C that co-orbital configurations that would be unstable on their own can be stabilised during the protoplanetary disc phase by the presence of another planet trapped in first order mean motion resonances either inside, or outside, of the co-orbital pair.

5. Evolution in a protoplanetary disc: 2D hydrodynamical simulations

The torques applied on the planets in the previous sections were obtained for a single planet embedded in a disc. When there is more than one planet in the disc, the surface density perturbations by the other planets can alter the torque on each individual planet (Baruteau & Papaloizou 2013; Pierens & Raymond 2014; Brož et al. 2018). On the other hand, moderate-mass planets might open a partial gap around their orbits that can also vary the torque from its pure type-I value. In this section we take these effects into account by running two-dimensional (2D) locally isothermal hydrodynamical simulations using FARGO¹ code (Masset 2000) for a system with two co-orbital planets in different mass regimes.

5.1. Disc and planets setups

The disc in our simulations is extended radially from 0.3 to 2.5 au and azimuthally over the whole 2π. It is gridded into N_r × N_ϕ = 873 × 1326 cells with logarithmic radial segments. The resolution is chosen such that the half horseshoe width of a 3 M_⊕ planet can be resolved by about 6 cells. The surface density profile is Σ = Σ₀r^−α = 2 × 10⁻⁴r^−0.85 in code units. This corresponds to 1777 g cm⁻² when the radial unit is 1 au and the mass unit is 1 M_⊙. The disc viscosity follows the alpha prescription of viscosity ν_visc = α_viscc_sH where c_s is the sound speed and H is the disc scale height. These two quantities are related to the aspect ratio h as h = H/r = c_s/v_k, v_k being the Keplerian velocity. The disc is flared with the aspect ratio of h = h₀r^f = 0.05r^0.175. The reason of such choices for surface density and aspect ratio profiles will be explained in Sect. 5.2.

Planets are initiated radially at r₁ = r₂ = 1 au and azimuthally at λ₁ = 0 (more massive one) and λ₂ = +50 or −50 degrees. The mass of the planets are increased gradually during the first 50 years to avoid abrupt perturbations in the disc while they are kept on circular orbits until t = 100 yr. The total time of the simulations is 2000 or 3000 years depending on the migration rate of the planets. In cases that the planets migrate faster, we had to stop the simulations earlier to avoid the effect of the inner boundary. In order to have the consistency with the 1D simulations, we used ϵ_p = 0.4h for smoothing the planet’s potential. The planets do not accrete gas during these simulations and their masses stay constant after t = 50 yr.

We used various planet mass pairs as given in the left columns of Table 1, where m₁ is the mass of the leading planet and m₂ the mass of the trailing one. We saved the forces and torque on each planet 20 times per year, allowing us to estimate the value of the different terms of the stability criterion u, given by Eq. (45). The method to compute these different terms is given in Appendix D. Columns 3–8 of Table 1 give the averaged values of each quantity over the whole simulation time, and the 9th column is the quantity 2u, the sum of these contributions.

5.2. Comparison of 1D to 2D discs

First, we compare the evolved disc in the hydro-models to our 1D discs. In Sect. 4.2, we saw that the criterion (56), based on the analytical torques from Tanaka et al. (2002), was a good approximation to estimate the stability of the system. This criterion is a function of the local flaring index and surface density slope f and α through the parameter K which parametrises the local slope of the torque felt by both planets. In all of the simulations presented in this section, the flaring index and initial surface density slopes are f = 0.175, α = 0.85 implying that the disc is viscously in equilibrium. These values correspond to K = 0.7, for which the co-orbital planets in type-I migration should always be diverging (see Fig. 7). In our hydrodynamical runs, we estimate K_j – the value of K which is felt by the jth planet – that can be linked to the quantities Γ_j0 and Γ_jΔ by the following relation (see Eq. (55)):

$$\begin{array}{c}\hfill K_j=(1+\frac{m_j{\mathrm{\Gamma }}_{j\mathrm{\Delta }}}{(m_1+m_2){\mathrm{\Gamma }}_{j0}})/2.\end{array}$$(60)

The K_js for the hydro simulations are hence computed and given in the last two columns of Table 1. It appears clearly that the local slope of the torque felt by each planet is different from what is expected from Tanaka et al. (2002). The upper block of the table (the first three rows) contains the models with low-mass planets such that they do not perturb the disc greatly. The maximum surface density perturbation δΣ/Σ₀ in these models is only about 3%. Therefore the difference of K_j with respect to the pure type-I migration which is expected for this type of planet is due to the presence of the co-orbital companion. For the models in the second block of the table, δΣ/Σ₀ is at most 35%. In these models, the torques are modified both by the presence of the other planet and the partial gap. Hence, the stability of co-orbitals is expected to be different from those obtained using a 1D disc model. However, we note that if one of the co-orbitals is significantly more massive than the other, the constant part of the torque that applies on it remains the same regardless of the position of the lower mass planet. For large mass discrepancies, the first two terms of u (Eq. (45)) might hence be properly estimated by the 1D model, see for example Fig. 11. However, hydrodynamical simulations are needed to estimates the Γ_jΔ and R_jζ terms.

5.3. Low to moderate mass planets: Super-Earth to mini-Neptune

In this section we study the evolution of co-orbitals in the super-Earth to mini-Neptune mass regime which are the planets that do not open a full gap in the disc. For comparison with the stability criterion developed in Sect. 3.3 (attractive equilibrium for negative u, repulsive for positive u), we give in Table 1 the quantity d(ζ_max − ζ_min)/dt, averaged over the simulation time. For the first two blocks of the table, where the no-gap and partial-gap opening planets are listed, the stability of the Lagrangian equilibrium is indeed correctly predicted by the criterion u as it has the same sign as d(ζ_max − ζ_min)/dt. On these sets of simulations, we can see a trend that was already remarked by Pierens & Raymond (2014): a more massive leading planet (m₁ >  m₂) tends to stabilise the co-orbital configuration, while if the more massive planet is trailing, the system slowly evolves away from the equilibrium.

To see how much the presence of the second planet can affect the torque, we present the torque analysis for the models (m₁, m₂) = (5, 3) M_⊕ in Fig. 13. In panel a, we plot the scaled torque Γ_j/Γ₀, with ${\mathrm{\Gamma }}_0={({\displaystyle \frac{m_j/m_0}{h}})}^2{\mathrm{\Sigma }}_{\mathrm{p}}{a}^4{\mathrm{\Omega }}_{\mathrm{p}}^2$, vs. distance from the planet’s orbit δr, which is scaled to the planet’s mutual Hill radius $R_{\mathrm{H}}=\sqrt[3]{(m_1+m_2)/3m_0}$. Γ_j/Γ₀ is the torque that is exerted on the jth planet by the material within ±δr of the planet’s orbit. In this figure, we compare the torque on each planet in the co-orbital simulation with the corresponding single-planet model. The scaled torques on single planet models are almost identical as expected for type-I migration. The torque on the co-orbital 5 M_⊕ follows the single planets up to about 5 R_H. It means that the torque from the co-rotation region and the its own spiral in this area is identical to the single planet model. As we move further out, the torque levels up until 8 R_H, decreases until 10 R_H, and increases again until it reaches the value of the total torque (dotted line). The cause of this variation can be found in panel b, where the surface density perturbation is shown, and in panel c, in which we plot the torque from each location in the disc on the planet. In panel c, the solid line demonstrates the torque from the inner disc and the dashed line is the negative of the torque from the outer disc. The advantage of this plot is that it shows where the torque from the disc in the co-orbital model differs from the single one. We see in panel c that the torque from the inner and outer disc is identical until about 5 R_H, where the positive torque from the inner disc increases in the co-orbital model. This area of the positive torque, marked by a red dot, corresponds to the area between the two inner arms of the planets. The presence of the second planet creates a slight depression in the surface density at the left side of the planet compared to the right. It makes the torque from this area more positive. As we move further the torque drops strongly due to the fact that we get close to the over-dense part of the second planet’s spiral arm which is located on the left side of the planet and exerts a negative torque on the 5 M_⊕. Then, this over-dense arm moves to the right side of the planet and its effect turns to a positive torque. The sum of all these components which arise from the depletion between the planets’ arms and other planet’s spiral is responsible for the deviation of the torque from the single planet models. The torque on the 3 M_⊕ is very similar except that the effect of the other planet is stronger due to its stronger arm (panel d): the diamond and square symbols mark the passage of the 5 M_⊕ planet’s spiral arms by the azimuthal position of the 3 M_⊕ planet.

As the spiral arms follow the planets in their libration around the Lagrangian equilibrium, both the torques and radial forces applied on each planet evolve over the libration time-scale, which is responsible for the non-negligible terms in Cols. 4, 5, 7, and 8 of Table 1.4pt

In Fig. 14, we present the same torque analysis for the model (m₁, m₂) = (6, 12) M_⊕, which has two partial gap opening planets. As in the low-mass co-orbital model, the main cause of the deviation from the single-planet models is the depression of the surface density between the two planets’ arms and the presence of the other planet’s arm. In addition, the partial gap between the two planets (ϕ ∈ [π, 4π/3]) is deeper than the rest of the gap, which creates an additional offset for the torques applied on each planet.

The torques on the co-orbital planets is hence very different from those that would apply on a single planet, and the difference originates from the suppression of the surface density between the planets and the effect of the other planet’s spiral arm. According to our current knowledge, there is no extensive study that tell us how the depletion of mass between the co-orbitals changes by the disc parameters or the mass of the planets. On the other hand, the torque from the other planet’s spiral arm depends on the strength of the arm which depends on the planet’s mass, and opening angle of the spiral which is a function of the disc aspect ratio. Hence, we expect the stability of the co-orbitals depends on these two parameters because as the planets librate around their equilibrium point, their distance from each other and consequently from each other’s spiral would also change. In the following section, we remove the complexity of the partial gap and the co-rotation torque by replacing one of the co-orbitals with a gap-opening planet.

5.4. Gap-opening planets: A Jupiter and an Earth

In order to isolate the effect of the spiral arm of one planet on the other, we ran two simulations with a Jupiter-mass planet and an Earth-mass planet. In one of the simulations the Earth-mass planet is leading, and in the other it is trailing.

In Fig. 15 we present the disc surface density perturbation and the torque analysis for these models. The torque on the Earth-mass planet in these models only comes from the Jupiter-mass planet, either from the material accumulated in its Hill radius or spirals. As the upper panel of Fig. 15 shows, the spirals of the Earth-mass planet are so weak because there is little material in the gap in the disc to form them. The lower panel shows that the sign of the torque on the Earth-mass planet depends on its location compared to the Jupiter-mass planet. Here we explain the torque analysis for the planet on the right side (indicated by a cross sign) and the opposite argument is applied for the planet on the left (marked by a plus sign). Following the solid blue line in the lower panel, we see the (negative) torque slowly increases until the red line, where the gap edge is located. This indicates that the main torque (see the dashed blue line) does not come from the material around the Jupiter-mass planet. As we add the contribution of the material from the red to the yellow line, a large negative torque is exerted on the planet by the inner spiral which is located to the left side of the planet. As we get further, the continuation of the inner spiral adds a positive torque but since it is farther than the section on the left, it cannot change the torque considerably. The grey line marked where we reached the inner edge of the disc, and therefore, the oscillations after the grey line only originate from the outer disc.

Based on the calculations in Sect. 3.3, the partial derivatives Γ_jΔ = ∂Γ_j/∂Δ and R_jζ = ∂R_j/∂(ζ − ζ_eq) are key parameters the for stability of the co-orbitals, see Cols. 4, 5, 7, and 8 of Table 1. Figure 16 represents the evolution over time of the torque and radial forces for the leading-Earth-mass planet case “×”. On the left panel, we can see that Γ_1Δ >  0, which leads to a positive term in the expression of u hence destabilising the configuration. On the other hand, the left panel shows that R_1ζ <  0, which induces a stabilising term. These two effects oppose one another and determine the attractiveness of the equilibrium. We note that in this particular case, our estimation of u did not match the evolution of the system: we computed a negative u, while the system is diverging from the equilibria. As the sum of the two dominant terms is an order of magnitude lower than each of these terms, this might be due to higher-order effects in the expansion of the torques and forces that were neglected when we computed the expression of u (in our runs, the semi-amplitude of libration is of ∼0.17 radians which is at the limit of the validity of the linear model).

The partial derivatives Γ_jΔ and R_jζ applied on the Earth-mass planet come from the spiral arm of the Jupiter-mass planet, see Fig. 15. As these effects are of opposite sign for leading and trailing planets (bottom panel of Fig. 15), they lead to a qualitatively similar behaviour for the stability around the L₄ and L₅ equilibria of the giant planet. To confirm this trend, we ran a set of simulations with a 10 M_⊕ planet either leading or trailing a Jupiter-mass planet for different disc parameters, varying the α parameter of the viscosity and the aspect ratio. Viscosity affects the gap depth and width, and the aspect ratio widens or tightens the spirals. For all the tested disc parameters, either both leading and trailing 10 M_⊕ planets converged toward the equilibria, or they both diverged away from it.

In Fig. 17 we show the results for these 7 different disc profiles that we tested. Seemingly, the stability of the planets inside the gap of a massive planet is a delicate trade-off between the disc parameters, although based on this small set, lower viscosity and smaller aspect ratios (that result in deeper and wider gaps) seem to stabilise co-orbital configurations. As this dependency is key to estimating the probability of the existence of Earth to super-Earth mass trojan companions to giant planets, we will investigate this topic more thoroughly in a future study.

6. Stability in the direction of the eccentricity and the inclinations

In this section we study the effect of dissipation on the evolution of the eccentricities and inclinations of the co-orbitals, for low values of e_j and I_j (≲0.1). At first order, the Poincaré variables (Eq. (1)) read:

$$\begin{array}{c}\hfill x_j=\frac{\sqrt{{\mathrm{\Lambda }}_j}}{\sqrt{2}}{e}_j{\mathrm{e}}^{i{\varpi }_j},\text{and}{y}_j=\frac{\sqrt{{\mathrm{\Lambda }}_j}}{\sqrt{2}}{I}_j{\mathrm{e}}^{i{\mathrm{\Omega }}_j}.\end{array}$$(61)

In the absence of dissipation, these variables follow the equations of variation given by the system (15).

We assume that the evolution of the orbital elements induced by dissipative forces can be modelled by migration and damping time-scales:

$$\begin{array}{c}\hfill {\dot{a}}_j=-a_j/{\tau }_{a_j},{\dot{e}}_j=-e_j/{\tau }_{e_j},\text{and}{\dot{I}}_j=-I_j/{\tau }_{I_j}.\end{array}$$(62)

We note that modelling the migration by such a law is equivalent to taking K = 0 in Sect. 4.1. It can be shown that the results of this section remain valid for any value of K, as the local variations of the torques over the resonant time-scale have a negligible effect on the evolution of the variable x_j and y_j. We hence consider the following non-conservative terms:

$$\begin{array}{cc}\hfill {\dot{x}}_{j,d}& =x_j(-\frac{1}{4{\tau }_{\mathit{aj}}}-\frac{1}{{\tau }_{\mathit{ej}}}),\hfill \\ \hfill {\dot{y}}_{j,d}& =y_j(-\frac{1}{4{\tau }_{\mathit{aj}}}-\frac{1}{{\tau }_{I,j}}),\hfill \end{array}$$(63)

for the evolution of the x_j and y_j. The equation of variations hence read:

$$\begin{array}{c}\hfill \dot{\mathit{x}}=M_x\mathit{x}+{\dot{\mathit{x}}}_d,\dot{\mathit{y}}=M_y\mathit{y}+{\dot{\mathit{y}}}_d,\end{array}$$(64)

where the M_x, M_y can be found in Appendix A, ${\dot{\mathit{x}}}_d=({\dot{x}}_{1,d},{\dot{x}}_{2,d})$ and ${\dot{\mathit{y}}}_d=({\dot{y}}_{1,d},{\dot{y}}_{2,d})$ (Eq. (63)). As we will be primarly interested in the evolution of the orbital elements e_j and I_j, we normalise the variables x_j: $X_j=e_j{\mathrm{e}}^{i{\varpi }_j}=x_j/\sqrt{{\mathrm{\Lambda }}_j/2}$ and $Y_j=I_j{\mathrm{e}}^{i{\mathrm{\Omega }}_j}=y_j/\sqrt{{\mathrm{\Lambda }}_j/2}$. The evolution of these new variables reads:

$$\begin{array}{cc}\hfill \dot{\mathit{X}}& =M_X(\zeta )\mathit{X}+{\dot{\mathit{X}}}_d,\hfill \\ \hfill \dot{\mathit{Y}}& =M_Y(\zeta )\mathit{X}+{\dot{\mathit{Y}}}_d,\hfill \end{array}$$(65)

and the dissipative terms read:

$$\begin{array}{cc}\hfill {\dot{X}}_{j,d}& =(\frac{{\dot{\mathrm{\Lambda }}}_j}{2{\mathrm{\Lambda }}_j}-\frac{1}{4{\tau }_{\mathit{aj}}}-\frac{1}{{\tau }_{\mathit{ej}}})X_j,\hfill \\ \hfill {\dot{Y}}_{j,d}& =(\frac{{\dot{\mathrm{\Lambda }}}_j}{2{\mathrm{\Lambda }}_j}-\frac{1}{4{\tau }_{\mathit{aj}}}-\frac{1}{{\tau }_{\mathit{Ij}}})Y_j.\hfill \end{array}$$(66)

6.1. Stability in the direction of the eccentricities

6.1.1. Constant masses

We study the stability in the direction x_j, related to the eccentricity and the argument of periastron, in the neighbourhood of the L₄ circular equilibrium for constant masses. In this section we consider dissipation time-scales that are not necessarily small with respect to $\frac{m_j}{m_0}{\eta }_{\mathrm{L}}$. The equation of variation of the variables x_j is given by

$$\begin{array}{c}\hfill \dot{\mathit{X}}=M_X(L_4)\mathit{X}+{\dot{\mathit{X}}}_d,\end{array}$$(67)

where M_X(L₄) is obtained by estimating the terms of M_X at the circular L₄ equilibria, given by Eq. (25). The system of Eq. (65) have two eigenvalues. At first order in ε:

$$\begin{array}{cc}\hfill g_{X\pm }& =i\epsilon \frac{g_{L_4}}{2}-\frac{\epsilon }{8}(\frac{m_1-m_2}{m_1+m_2}\frac{1}{{\tau }_{a-}}+\frac{4}{{\tau }_{e+}})\hfill \\ \hfill & \pm \epsilon \sqrt{{\left(\frac{1}{8{\tau }_{X-}}\right)}^2-i\frac{1}{2}\frac{m_1-m_2}{m_1+m_2}\frac{g_{L_4}}{{\tau }_{X-}}+{\left(i\frac{g_{L_4}}{2}\right)}^2},\hfill \end{array}$$(68)

where

$$\begin{array}{c}\hfill 1/{\tau }_{a\pm }=\frac{{\tau }_{a_1}\pm {\tau }_{a_2}}{{\tau }_{a_1}{\tau }_{a_2}}1/{\tau }_{e\pm }=\frac{{\tau }_{e_1}\pm {\tau }_{e_2}}{{\tau }_{e_1}{\tau }_{e_2}},\end{array}$$(69)

1/τ_X− = 1/τ_a− + 4/τ_e−, and

$$\begin{array}{c}\hfill g_{L_4}=\epsilon (-\frac{27}{8}+\frac{87}{16}\sqrt{3}{z}_{L_4})\frac{m_1+m_2}{m_0}{\eta }_{\mathrm{\Gamma }}.\end{array}$$(70)

In the conservative case, 1/τ_e± = 1/τ_a± = z_L4 = 0, and we obtain g₋ = 0 and g₊ = i27/8(m₁ + m₂)/m₀η. The direction associated to these eigenvalues were described in Sect. 2.3.1: the eccentric Lagrangian equilibrium, where e₁ = e₂ and ϖ₁ − ϖ₂ = ζ = ±π/3; and the anti-Lagrangian equilibrium, where m₁e₁ = m₂e₂ and ϖ₁ − ϖ₂ = ζ + π = ∓2π/3.

Figure 18 shows the values of τ_e2/τ_a1 for which the real component of the eigenvalues (68) vanishes, with respect to τ_a1/τ_a2. These plots were made using m₁ = 10⁻⁴m₀, τ_e1 = τ_e2m₂/m₁, and τ_a1 = 10/m₁. For a given mass ratio, the manifold e₁ = e₂ = 0 is attractive below the two curves of the given colour. Above the solid line, the system diverges following the anti-Lagrangian direction; while systems above the dashed curve diverge following the eccentric Lagrangian direction. These curves were obtained using several assumptions on the relations between the masses and the damping and migration time-scales, and that the stability of the X_j directions in the τ_a1/τ_a2, τ_e2/τ_a1 plane depends greatly on these assumptions.

Figure 19 represents the evolution of two configurations taken on the left border of Fig. 18: τ_a1/τ_a2 = 0.01. In the left panel, τ_e2/τ_a1 = 5, while on the right τ_e2/τ_a1 = 20. The stability in the X_j directions is given by the position of the configurations relative to the blue curves of Fig. 18. In both cases, the motion relative to the direction of the eccentric Lagrangian equilibria (black solid lines in Fig. 19) is quickly damped as we are far below the dashed lines in both cases. As the quantity ϖ₁ − ϖ₂ converges toward 240 ° ( = ζ + 180°), which is the direction of the anti-Lagrangian equilibria, the eccentricity either decreases as this direction is stable (the left case is below the solid blue curve of Fig. 18), or increases if the anti-Lagrangian direction is unstable (the right case is above the solid blue curve).

6.1.2. Effect of mass change

We now consider the effect of a slow, isotropic mass change on the eccentric Lagrangian equilibria and anti Lagrangian equilibria previously discussed. For this mass change to impact the evolution of the configuration, it has to be comparable to the migration and damping time-scales. We hence assume that the perturbative terms 1/τ_aj, 1/τ_ej and ṁ_j/m_j are of size ε². The details of the computations can be found in Appendix E.

• Evolution of the eccentricities along the anti-Lagrangian equilibria: taking 𝒳₂ = 0, Eq. (E.5) yields, at second order in ε:

  $$\begin{array}{c}\hfill {\dot{e}}_j={\epsilon }^2\frac{m_k{\mathcal{X}}_1{\overline{\mathcal{X}}}_1}{4e_j(m_1+m_2)}(4m_j{\dot{m}}_k-4m_k{\dot{m}}_j+\frac{m_k}{T_{AL4}})\end{array}$$(71)

  where ${\mathcal{X}}_1{\overline{\mathcal{X}}}_1$ is a positive real quantity and

  $$\begin{array}{c}\hfill \frac{1}{T_{AL4}}=m_2(-\frac{1}{{\tau }_{a1}}+\frac{1}{{\tau }_{a2}}-\frac{4}{{\tau }_{e1}})+m_1(-\frac{1}{{\tau }_{a2}}+\frac{1}{{\tau }_{a1}}-\frac{4}{{\tau }_{e2}})\end{array}$$(72)

• Evolution of the eccentricities along the eccentric-Lagrangian equilibria: taking 𝒳₁ = 0, Eq. (E.5) yields, at second order in ε:

  $$\begin{array}{c}\hfill {\dot{e}}_j=-{\epsilon }^2\frac{{\mathcal{X}}_2{\overline{\mathcal{X}}}_2}{e_j(m_1+m_2)}(\frac{m_1}{{\tau }_{e1}}+\frac{m_2}{{\tau }_{e2}})\end{array}$$(73)

  where ${\mathcal{X}}_2{\overline{\mathcal{X}}}_2$ is a positive real quantity.

If the eccentricities are damped by the disc (τ_ej >  0) then the mode associated to the Eccentric Lagrangian equilibrium will always be damped toward e₁ = e₂ = 0. However the eccentricities can increase along the Anti-Lagrangian equilibria if the more massive of the two planets migrate inward (or accrete gas) fast enough (see Eq. (71)).

6.2. Stability in the direction of the inclinations

6.2.1. Constant masses

We study the stability in the direction Y_j, related to the inclinations and the ascending nodes of the co-orbitals, for any amplitude of libration of the resonant angle. Eq. (A.2) becomes:

$$\begin{array}{c}\hfill M_Y=\left(\begin{array}{cc}\frac{{\overline{A}}_Y}{m_1}-Y_{1,d}/Y_1& \frac{-{\overline{A}}_Y}{\sqrt{m_1{m}_2}}\\ \frac{-{\overline{A}}_Y}{\sqrt{m_1{m}_2}}& \frac{{\overline{A}}_Y}{m_2}-{\dot{Y}}_{2,d}/Y_2\end{array}\right).\end{array}$$(74)

At first order in ε, M_Y(L₄) can be diagonalised, with the diagonal elements being:

$$\begin{array}{cc}\hfill s_{Y\pm }& =i\frac{g_Y}{2}-\frac{1}{8}(\frac{m_1-m_2}{m_1+m_2}\frac{1}{{\tau }_{a-}}+\frac{4}{{\tau }_{I+}})\hfill \\ \hfill & \pm \sqrt{{\left(\frac{1}{8{\tau }_{Y-}}\right)}^2-i\frac{1}{2}\frac{m_1-m_2}{m_1+m_2}\frac{g_{L_4}}{{\tau }_{Y-}}+{\left(i\frac{g_{L_4}}{2}\right)}^2},\hfill \end{array}$$(75)

where

$$\begin{array}{c}\hfill 1/{\tau }_{a\pm }=\frac{{\tau }_{a_1}\pm {\tau }_{a_2}}{{\tau }_{a_1}{\tau }_{a_2}}1/{\tau }_{I\pm }=\frac{{\tau }_{I_1}\pm {\tau }_{I_2}}{{\tau }_{I_1}{\tau }_{I_2}},\end{array}$$(76)

1/τ_X− = 1/τ_a− + 4/τ_I−, and

$$\begin{array}{c}\hfill g_Y=\frac{{\overline{A}}_Y}{2}\frac{m_1+m_2}{m_0}{\eta }_{\mathrm{\Gamma }}.\end{array}$$(77)

We note that the eigenvalues (75) have a similar expression to that in the direction of the eccentricity, Eq. (68), but here the results are valid for any amplitude of libration in the trojan and horseshoe domains. The amplitude of libration affects the stability in the Y_j direction through the value of ${\overline{A}}_Y$ (see Fig. 3).

The orbits can either be attracted toward I₁ = I₂ = 0, or diverge following m₁I₁ = m₂I₂ and Ω₂ = Ω₁ + π, or I₂ = I₁ and Ω₂ = Ω₁, see Sect. 2.3.2. Examples of convergence and divergence along m₁I₁ = m₂I₂, Ω₂ = Ω₁ + π are shown in Fig. 20.

6.2.2. Effect of mass change

We now consider the effect of a slow, isotropic mass change on the inclination of quasi-circular co-orbitals (up to first order in eccentricities). As we did for the study of the evolution of the eccentricities, we assume that the perturbative terms 1/τ_aj, 1/τ_Ij and ṁ_j/m_j are of size ε², and that the mass evolution is isotropic. The details of the computations are identical to the eccentric case and can be found in Appendix E.

• Evolution of the inclinations along m₁I₁ = m₂I₂, Ω₂ = Ω₁ + π: taking 𝒴₂ = 0, we obtain, at second order in ε:

  $$\begin{array}{c}\hfill {\dot{I}}_j={\epsilon }^2\frac{m_k{\mathcal{Y}}_1{\overline{\mathcal{Y}}}_1}{4I_j(m_1+m_2)}(4m_j{\dot{m}}_k-4m_k{\dot{m}}_j+\frac{m_k}{T_I})\end{array}$$(78)

  where ${\mathcal{Y}}_1{\overline{\mathcal{Y}}}_1$ is a positive real quantity and

  $$\begin{array}{c}\hfill \frac{1}{T_I}=m_2(-\frac{1}{{\tau }_{a1}}+\frac{1}{{\tau }_{a2}}-\frac{4}{{\tau }_{I1}})+m_1(-\frac{1}{{\tau }_{a2}}+\frac{1}{{\tau }_{a1}}-\frac{4}{{\tau }_{I2}})\end{array}$$(79)

• Evolution of the inclination along I₁ = I₂, Ω₁ = Ω₂: taking 𝒴₁ = 0, we obtain, at second order in ε:

  $$\begin{array}{c}\hfill {\dot{I}}_j=-{\epsilon }^2\frac{{\mathcal{Y}}_2{\overline{\mathcal{Y}}}_2}{I_j(m_1+m_2)}(\frac{m_1}{{\tau }_{I1}}+\frac{m_2}{{\tau }_{I2}})\end{array}$$(80)

  where ${\mathcal{Y}}_2{\overline{\mathcal{Y}}}_2$ is a positive real quantity.

If the inclination are damped by the disc (τ_Ij >  0) then the mode associated to I₁ = I₂, Ω₁ = Ω₂ will always be damped toward I₁ = I₂ = 0. However the inclinations can increase along m₁I₁ = m₂I₂, Ω₂ = Ω₁ + π if the more massive of the two planets migrate inward (or accrete gas) fast enough (see Eq. (78)). We remind the reader that these results on the inclinations are valid for any amplitude of libration of the resonant angle, up to horseshoe orbits.

7. Summary and conclusions

7.1. Summary

In this paper we have studied the stability of the co-orbital resonance under dissipation in the planetary case ((m₁, m₂)≪m₀). In Sect. 3 we developed an integrable model of the 1:1 MMR perturbed by a generic dissipation and derived the stability conditions of the L₄ and L₅ equilibria.

In Sect. 3.2 we showed that under the effect of a constant torque applied on each planet, the phase space of the resonance becomes asymmetric, as the position of the Lagrangian equilibria L₃, L₄ and L₅ change. The tadpole (trojan) area is larger if the torque per mass unit applied on the leading planet is greater than the torque per mass unit applied on the trailing one. We also saw that if the difference between these two torques is too large, two out of the three equilibrium points could merge and vanish, leading to a phase space with a single equilibrium point. These results are in agreement with those of Sicardy & Dubois (2003), obtained in the restricted case (m₁ ≪ m₀,  m₂ = 0). This effect can also contribute to the instability observed by Pierens & Raymond (2014), where they showed that similar mass co-orbitals were unstable during the partial gap-opening regime, due to the opposite torques induced by a higher gas depletion between the two planet than everywhere else in the gap.

In Sects. 3.2.2 and 3.3, we then studied the stability of the Lagrangian equilibria L₄ and L₅ as a function of the forces applied on each planet, their masses m_j, and the evolution of their mass ṁ_j. This study can be split into two parts:

• First, the evolution of the masses of the co-orbitals, along with the constant torques that are applied on them, change the width of the co-orbital resonances. It can lead to either a convergence toward the Lagrangian equilibria in the case of outward migration (positive total torque, Γ₁₀ + Γ₂₀ >  0) or overall mass increase (ṁ₁ + ṁ₂ > 0), while inward migration (negative total torque) and mass loss induces a slow divergence from the Lagrangian equilibria. These results are in agreement with those of Fleming & Hamilton (2000), which were obtained in the restricted case (m₁ ≪ m₀,  m₂ = 0).

• Second, if the forces applied on each planet vary over the resonant time-scale, we show that the dependency of the torques on the semi-major axis, and the dependency of the radial component of the perturbative forces on the value of the resonant angle, impact significantly the stability of the system (the effect of a radial dependency of the torque was discussed by Sicardy & Dubois 2003, in the restricted case).

These two effects were considered to derive the stability criterion u for the Lagrangian equilibria (Eq. (45)).

Sections 4 and 5 were dedicated to comparing these results to N-body simulations in 1D disc models, and hydrodynamic simulations. In Sect. 4, we applied type-I migration prescriptions on a pair of planets in an evolving protoplanetary disc. The stability criterion successfully predicts the stability of the system, as a function of their masses, their migration time-scale τ_aj, and the slope of that migration parametrised by K_j. In addition, running planetary system evolution through the disc lifetime allowed us to study the balance between the destabilising effect of inward migration and the stabilising effect of mass accretion: first, planets tend to grow in mass significantly faster than they migrate, which leads to a convergence toward the exact equilibrium. However, in the later stages of the disc lifetime, the planets migrate quickly, leading to a divergence from the equilibrium. In addition, we showed that co-orbitals that belong to a resonant chain with other planets can be stabilised during the migration phase.

However, the comparison to hydrodynamics simulations show the limits of the 1D models: despite having similar initial conditions for the disc, the forces that are applied on each planet in the hydrodynamical simulation are totally different from those given by type-I prescriptions. Indeed, as the two planets evolve around the same semi-major axis, the disc is significantly perturbed both radially and azimuthally (Fig. 13, see also Brož et al. 2018). It creates structures whose effects cannot be azimuthally averaged, as they follow the position of the planets. Notably, as both planets librate around the Lagrangian equilibria, they move relatively to one another’s spiral arms. The additional torques and radial forces applied on each planet hence evolve over the libration time-scale, that can either have a stabilising effect, or destabilising one, see Table 1. It is the sum of all these terms that dictates the evolution of the system.

In the super-Earth range (3–5 M_⊕) we note a trend that was observed by Pierens & Raymond (2014): more massive leading planets tend to stabilise the system. We show here that this stabilisation is due to the variations of the torques felt by each planet over the resonant libration, as they are successively closer to, and then further away from, one another’s spiral arm. This trend is also present in the mini-Neptune regime (up to 15 M_⊕) with the apparition of other structures, such as a partial gap that is deeper between the co-orbitals.

Finally, in the case where the gap is totally open, we ran a set of simulations with a Jupiter-mass planet trailed or preceded by Earth or super-Earth mass planets. Here the dominant effect for the stability was the variation of the radial forces and torques applied on the Earth-mass planet by the Jupiter-mass planet’s spiral arms, during each libration period. The symmetry of the spiral arms with respect to the Jupiter-mass planet led to a similar behaviour for leading and trailing smaller mass companions: for all tested disc profiles, both leading and trailing companions behaved in a similar way (both diverging from or both converging toward L₄/L₅). However, as shown in Fig. 17, different disc parameters change the stability of such configurations. The effect of the disc parameters on the shape and strength of the Jupiter-mass planet’s spiral arms will be the subject of a future study.

In Sect. 6, we studied the stability of the Lagrangian equilibrium in the direction of the eccentricities (at first order), and the stability of the whole tadpole and horseshoe domain in the direction of the inclinations (at first order as well). We have shown that even in the case were the dissipative forces tend to damp the eccentricities and inclination of the planets, those could increase along a particular family of orbits.

7.2. Conclusions

7.2.1. On the limitation of 1D disc models

We have shown that disc-planet coupling generates structures in the disc that cannot be azimuthally averaged, leading to variations over time of the torques and radial forces that applies on each planet. Using similar disc profiles in 1D and hydrodynamical simulations, these differences lead to opposite results on the stability of the Lagrangian equilibria. Similar observations were made by Brož et al. (2018) in a more general context.

7.2.2. On the evolution of co-orbitals

Trojan swarms. The asymmetry between the L₄ and L₅ domains induced by the difference of torque per mass unit (Fig. 4) can be used to explain the potential asymmetry between the leading and trailing Jupiter’s and Neptune’s trojan swarm However it requires one to properly estimate the torques that are felt by each of the asteroids: we showed in Sect. 5 that the torque per mass unit applied by the protoplanetary disc on the leading and trailing Trojans are not negligible, and comes mainly from the Jupiter-mass planet’s spiral arms (see Fig. 15). These torques will hence strongly depend on the disc parameters, and are of opposite sign for L₄ and L₅ Trojans. As a result, L₄ and L5 domains would be more symmetric than if we apply the same torque on all asteroids.

Co-orbital exoplanets. We have shown that the attractiveness of the Lagrangian equilibria depends on the mass distribution between the planet, the total mass, the accretion rate, the constant torques and radial forces that apply on each planet, but also on how these quantities evolve on the resonant time-scale. We have shown that long inward migration destabilises co-orbitals, while outward migration and mass accretion tend to stabilise them. Figure 11 shows that the stabilising terms coming from mass accretion is comparable to the destabilising terms coming from inward migration, and hence both have to be taken into account to properly estimate the stability of a system. However, this stabilising effect comes into play mainly in the earlier phase of the planet’s evolution. While in the later stages, its evolution is dominated by the migration.

As in Pierens & Raymond (2014), we also found that leading massive trojans tend to stabilise the configuration. In their study, these authors also showed that equal mass co-orbitals can be disrupted during the gap opening stages. We have shown that the stability of Earth-mass planets as trojan companions of a Jupiter-mass planet depend on the disc parameters, but that both L₄ and L₅ configurations tend to be stable or unstable for a given set of disc parameters.

We have shown that unstable co-orbital configurations could be stabilised by being trapped in first order mean motion resonance with a third planet, although in this part of the study we neglected the perturbation coming from the different planet’s spiral arms (Brož et al. 2018).

It is worth noting that the Lagrangian equilibria being repulsive does not necessarily imply that no co-orbital configurations can remain, it only implies that the amplitude of libration around the Lagrangian equilibria slowly increases over the migration time-scale, although that can lead to trojan orbits becoming horseshoe orbits, or even exiting the resonance. Similarly, attractive Lagrangian equilibria only implies a slow convergence toward it, but the configuration can still be disrupted on shorter time-scales for example through N-body interaction with other planets (Robutel & Bodossian 2009; Leleu et al. 2019).

7.2.3. On the detectability of co-orbitals exoplanets

In our hydro-simulation runs, in the [3, 15] M_⊕ range and for a given disc profile, all configurations with a leading more massive planet were attracted toward the Lagrangian equilibria for planets. On the contrary, for the Jupiter-mass planet’s Earth sized trojan the stability seemed to depend very little on who is leading in the orbit, but we showed that different disc parameters can change the attractiveness of the Lagrangian equilibria. In addition, Cresswell & Nelson (2009) found that during the co-orbital’s evolution in the disc, the mass discrepancy between the two planets keep increasing because the more massive planet starves-off the other.

Our study of the stability of the Lagrangian equilibria in the inclined direction also leads to important conclusions regarding the detectability of co-orbitals. We have shown that as long as the disc tend to damp inclinations, the system can evolve toward two directions: either coplanar co-orbitals, or mutually inclined co-orbitals following the m₁I₁ = m₂I₂, Ω₁ = Ω₂ + π direction. This later direction is favoured if the proper migration of the more massive of the two planets, or its mass accretion rate, is faster than the inclination damping of the smaller planet. As the inclination damping of the smaller planet is reduced by the deeper partial or full gap created by the more massive planet, that could significantly reduce the transit probability of both co-orbitals. Similarly, even when the disc damps the eccentricities of the two planets, these eccentricities can increase following the anti-Lagrangian equilibria m₁e₁ = m₂e₂ω₁ − ω₂ = ζ + π (Giuppone et al. 2010; Leleu et al. 2018).

Mutually inclined co-orbitals can still be detected using transit timing variations (TTVs; Ford & Holman 2007; Vokrouhlický & Nesvorný 2014; Leleu et al. 2019) or radial velocities (Laughlin & Chambers 2002; Leleu et al. 2015), however, these methods require that the co-orbitals librate with a significant amplitude around the Lagrangian equilibrium, and that the observations baseline is at least comparable with the libration time-scale. In addition, the planets have to be of comparable masses for the radial velocity method, as well as for TTVs if it is the larger of the two planets that is transiting. Finally, even in the absence of libration, the combination of transit and radial velocity measurements can be used to detect co-orbital configurations (Ford & Gaudi 2006; Leleu et al. 2017), although this requires good constraints on the eccentricity of the transiting planet.

---

Acknowledgments

The authors acknowledge support from the Swiss NCCR PlanetS and the Swiss National Science Foundation. S.Ataiee acknowledges the support of the DFG priority program SPP 1992 “Exploring the Diversity of Extrasolar Planets (KL 650/27-1)”.

References

Appendix A: Equations of the coorbital resonance at first order in e and I

The averaged Hamiltonian of the circular coplanar coorbital resonance is (Robutel & Pousse 2013):

$$\begin{array}{cc}\hfill {\overline{\mathcal{H}}}_0& ={\mu }_0(\frac{m_1^3}{Z^2}+\frac{m_2^3}{{(Z_2-Z)}^2}\hfill \\ \hfill & +\mathcal{G}{m}_1{m}_2[\frac{m_1{m}_2}{{\mathrm{\Lambda }}_1{\mathrm{\Lambda }}_2}cos\zeta -(\frac{Z^4}{m_1^4}\hfill \\ \hfill & {+\frac{{(Z_2-Z)}^4}{m_2^4}-2\frac{Z^2{(Z_2-Z)}^2}{m_1^2{m}_2^2}cos\zeta )}^{-1/2}]).\hfill \end{array}$$(A.1)

While the equation of variations of the x and y variables are given by (Robutel & Pousse 2013; Robutel et al. 2015):

$$\begin{array}{c}\hfill \dot{\mathit{x}}=M_x(\zeta )\mathit{x},\dot{\mathit{y}}=M_y(\zeta )\mathit{y},\end{array}$$(A.2)

with

$$\begin{array}{c}\hfill M_x(\zeta )=\left(\begin{array}{cc}\frac{A_x(\zeta )}{m_1}& \frac{{\overline{B}}_x(\zeta )}{\sqrt{m_1{m}_2}}\\ \frac{B_x(\zeta )}{\sqrt{m_1{m}_2}}& \frac{A_x(\zeta )}{m_2}\end{array}\right),M_y=\left(\begin{array}{cc}\frac{A_y(\zeta )}{m_1}& \frac{{\overline{B}}_y(\zeta )}{\sqrt{m_1{m}_2}}\\ \frac{B_y(\zeta )}{\sqrt{m_1{m}_2}}& \frac{A_y(\zeta )}{m_2}\end{array}\right),\end{array}$$(A.3)

with

$$\begin{array}{cc}\hfill A^{(v)}& =-i\frac{m_1{m}_2}{2m_0}\eta (1-\frac{1}{\delta {(\zeta )}^3})cos\zeta ,\hfill \\ \hfill B^{(v)}& =i\frac{m_1{m}_2}{2m_0}\eta (1-\frac{1}{\delta {(\zeta )}^3}){exp}^{i\zeta },\hfill \\ \hfill A^{(h)}& =\frac{1}{4\delta {(\zeta )}^5}(5cos2\zeta -13+8cos\zeta )-cos\zeta ,\hfill \\ \hfill B^{(h)}& ={exp}^{-2i\zeta }-\frac{1}{8\delta {(\zeta )}^5}({exp}^{-3i\zeta }+16{exp}^{-2i\zeta }\hfill \\ \hfill & -26{exp}^{-i\zeta }+9{exp}^{i\zeta }).\hfill \end{array}$$(A.4)

Appendix B: Stability of partial equilibria

We can study the stability of the Lagrangian points even if the equations of variation (38) are not constant over time by studying the stability of partial equilibria (Vorotnikov 2002). To do so, we divide the variables two groups: the variables with respect to which the stability is investigated z = (z₁, z₂), and the remaining variable γ = (L, m₁,m₂). Their equations of variation is given by the system:

$$\begin{array}{c}\hfill \{\begin{array}{c}\frac{{\dot{z}}_j}{z_j}=u(L,m_1,m_2)+{(-1)}^j\nu (L,m_1,m_2),\hfill \\ \dot{L}={\mathrm{\Gamma }}_1+{\mathrm{\Gamma }}_2+\frac{{\dot{m}}_1+{\dot{m}}_2}{m_1+m_2}L,\hfill \\ \frac{m_j}{m_j}={\mathrm{Acc}}_j(L,m_1,m_2).\hfill \end{array}\end{array}$$(B.1)

where Acc_j is the accretion rate of the planet j. The system (B.1) can be rewritten:

$$\begin{array}{c}\hfill \dot{\mathit{z}}=\mathit{Z}(t,\mathit{z},\mathit{\gamma }),\dot{\mathit{\gamma }}=F(t,\mathit{z},\mathit{\gamma }),\mathit{Z}(t,\mathbf{0},\mathit{\gamma })=0.\end{array}$$(B.2)

We note that the equations of variations of each component of the vector z are uncoupled (Eq. (B.1)). We hence study the stability in the direction of each component z_j separately.

Following (Vorotnikov 2002), let a(r) and b(r) be arbitrary continuous, monotone increasing functions for r ∈ [0, h], where h is a positive real number, and such as a(0) = b(0) = 0. If for the system (B.2) a scalar function V exist such that

$$\begin{array}{c}\hfill a(||{\mathit{z}}_j||)\le V(t,{\mathit{z}}_j,\gamma )\le b(||{\mathit{z}}_j||)(a),\\ \hfill \dot{V}\le 0(b)\text{and}V(t,\mathbf{0},\mathbf{0})\equiv 0(c),\end{array}$$(B.3)

then the set z_j = 0 is uniformly stable. To verify these conditions, we simply take:

$$\begin{array}{c}\hfill V(t,{\mathit{z}}_j,\mathrm{\Gamma })=a(||{\mathit{z}}_j||)=b(||{\mathit{z}}_j||)=||{\mathit{z}}_j{||}^2.\end{array}$$(B.4)

(a) and (c) are automatically verified, and as $\dot{V}=2Re(u\pm i\nu ){\mathit{z}}_j{\overline{\mathit{z}}}_j$, (b) is verified if u is negative or null. The Lagrangian point L₄ is hence uniformly stable if u is negative or null.

Appendix C: Effect of resonant chains

The existence of trojans with small amplitudes of libration in the synthetic planetary systems described in Sect. 4.3 prompted us to study the effect of resonant chains on the co-orbital configuration. A complete study of the stability of co-orbitals in resonant chains is beyond the scope of this paper, so we restrain this analysis to the effect of a single planet inside or outside of a co-orbital configuration, trapped in a 3:2 or 4:3 mean motion resonance (MMR).

To begin with, we look at three examples where two co-orbitals (m₁/m₂ = 6) are captured in a 4:3 MMR with an outer planet m₃ such that m₃/m₁ = [1, 1.6, 2.5]. Figure C.1 shows the evolution of the amplitude of libration as well as the eccentricity of the co-orbitals after their capture. In these three cases, the amplitude of libration, which was initially at ζ_max − ζ_min = 320°, is quickly reduced down to a few tens of degree by the capture into the 4:3 MMR with the 3rd planet. However, past the first 2 × 10⁵ years, the chosen examples exhibit three different behaviours: the blue case sees its amplitude of libration monotonically decrease over time, while the purple one keeps increasing after the first phase of the capture. In both of these cases, the eccentricities of the co-orbitals reach an equilibrium value, typical for two planets migrating in a first order MMR. The red case displays a more complex behaviour for both its amplitude of libration and eccentricities.

The effect of the relative masses between the co-orbitals m₁ and m₂, and the 3rd planet m₃ is studied by integrating a grid of initial conditions, taking for the co-orbitals a₁ = a₂ = 1 au, ζ = 20° (which yield an amplidue of libration of 320°, hence a horseshoe configuration), and masses m₁/m₂ in the [1:20] range. Both the 3:2 and 4:3 MMR are studied. The mass of the 3rd planet is m₃ = 3 × 10⁻⁵m₀ when it is at a larger semi-major axis than the co-orbitals, and m₃ = 1 × 10⁻⁵m₀ when it is at a smaller semi-major axis. All eccentricities and inclinations are initially set to 0. τ_aj = 10/m_j, and τ_ej = τ_aj/150, which are consistent with the parameters of the disc described in Sect. 4.2 and Tanaka & Ward (2004). We set K = 0 (τ_aj independent from a_j). Each initial condition (each set of masses) is integrated for 1.6 × 10⁵ years, which corresponds to a migration down to ≈0.4 au for the co-orbitals, depending on the chosen masses.

Results are displayed in Figs. C.2 and C.3. In both figures white pixels represent the systems for which the co-orbitals were not in the intended resonance with m₃ at the end of the simulation. For the systems close to m₃ = m₁, the planet did not converge or did not converge fast enough to reach the desired resonance during the simulation. For the other white pixels, the desired resonance was crossed but the capture did not happen, or did not hold.

In Fig. C.2, each pixel represents the final amplitude of libration of the co-orbital configuration for that set of masses. In almost all the studied cases, the capture into a mean-motion resonance with another planet led the co-orbital configuration to greatly reduce its own amplitude of libration, going from horseshoe to trojan configuration.

Figure C.3 shows the evolution of the amplitude of libration of the co-orbitals once they are captured in the MMR with m₃. The set of integrations is the same as in Fig. C.2. The quantity −d(ζ_max − ζ_min)/da is obtained by comparing the amplitude of libration of the co-orbitals between t = [t_max/3 : 2t_max/3] and [2t_max/3 : t_max]. Areas of blue or purple pixels follow a similar behaviour to the blue example of Fig. C.1, where the amplitude of libration keeps decreasing post-capture. Orange or yellow pixels show mass ratios for which the amplitudes of libration are, on average, increasing post-capture. For comparison, in absence of m₃ these co-orbital configurations would be diverging from the equilibrium for all values of m₁/m₂ = (1 : 20], since in these examples, K = 0 (see Fig. 7). The dependency of the stability with respect to the mass ratios is obviously complex and will be the object of a future study. We can nonetheless see from Figs. C.2 and C.3 that resonant chains can have a stabilising effect on the co-orbital configuration.

Appendix D: Forces partial derivatives

The forces applied by the disc on each planet during the hydrodynamical simulations were saved with a time-step of 0.05 year. To compute the partial derivatives ${\mathrm{\Gamma }}_{j\mathrm{\Delta }}=\frac{\partial {\mathrm{\Gamma }}_j}{\partial \mathrm{\Delta }}$ and $R_{j\zeta }=\frac{\partial R_j}{\partial \zeta }$ we first performed a sliding averaging of the quantities Γ_j and R_j over the local orbital period, in order to get rid of short terms effect, notably the oscillations due to the eccentricities. The raw torque and the result of this averaging are displayed in black and red in Fig. D.1, respectively, in the case m₁ = 6 M_⊕ (leading), m₂ = 12 M_⊕ (trailing), discussed in Sect. 4.

Then we computed Γ_j0 and R_j0 by performing an additional sliding average over the local libration period. This period was computed by frequency analysis of the resonant angle ζ. For our example, the instantaneous value of Γ₂₀ is displayed in blue in Fig. D.1. The quantity Γ_j − Γ_j0 is shown in Fig. D.2, with respect to the variable Δ and time (colour code). As we made a linear approximation in the neighbourhood of Δ = z = 0 in the analytical part of the study, we fitted Γ_jIΔ + Γ_jζz to the quantity Γ_j − Γ_j0. Here again, this fit is done over a sliding window of width 1-libration period.

Appendix E: Evolution of eccentricity under mass change

In the conservative case, described in Sect. 2.3.1, the diagonalisation is obtained by a change of variable X = P_X𝒳, where the columns of the matrix P_X are proportional to the eigenvectors of the matrix M_X (this change is thus not unique). In the case of evolving masses, the matrix P_X is not constant. We obtain the following relation:

$$\begin{array}{c}\hfill \dot{\mathcal{X}}=M_{\mathcal{X}}\mathcal{X}=(P_X^{-1}{M}_X{P}_X-P_X^{-1}{\dot{P}}_X)\mathcal{X}.\end{array}$$(E.1)

We hence look for a change of variables 𝒳 = P_𝒳X that diagonalises M_𝒳 at second order in ε. To do so, we look for a change of basis ε close to P_X: $P_{\mathcal{X}}=P_X+\epsilon P{\prime }_X$. Noting P_𝒳[j, k] the kth element of the jth line, the chosen change of basis is:

$$\begin{array}{cc}& P_{\mathcal{X}}[1,1]=m_2{\mathrm{e}}^{i\frac{\pi }{3}}\hfill \\ \hfill & P_{\mathcal{X}}[1,2]={\mathrm{e}}^{i\frac{\pi }{3}}+\frac{\epsilon {\mathrm{e}}^{-i\frac{\pi }{6}}}{g_{L_4}}(\frac{{\dot{X}}_{1,d}}{X_1}-\frac{{\dot{X}}_{2,d}}{X_2})\hfill \\ \hfill & P_{\mathcal{X}}[2,1]=-m_1+\frac{\epsilon i}{g_{L_4}}[{\dot{m}}_1-\frac{m_1}{m_2}{\dot{m}}_2+m_1(\frac{{\dot{X}}_{1,d}}{X_1}-\frac{{\dot{X}}_{2,d}}{X_2})]\hfill \\ \hfill & P_{\mathcal{X}}[2,2]=1\hfill \end{array}$$(E.2)

while M_𝒳 reads:

$$\begin{array}{c}\hfill {\dot{M}}_{\mathcal{X}}=\left(\begin{array}{cc}{g}_1+r_{x1}& 0\\ 0& g_2+r_{x2}\end{array}\right)\left(\begin{array}{c}{z}_1\\ z_2\end{array}\right)\end{array}$$(E.3)

where

$$\begin{array}{cc}& g_1=\epsilon i\frac{27(m_1+m_2){\eta }_{\mathrm{L}}}{8m_0},\hfill \\ \hfill & r_{x1}={\epsilon }^2\frac{m_2\frac{{\dot{X}}_{1,d}}{X_1}+m_1\frac{{\dot{X}}_{2,d}}{X_2}-{\dot{m}}_1-{\dot{m}}_2}{m_1+m_2},\hfill \\ \hfill & g_2=0,\hfill \\ \hfill & r_{x2}={\epsilon }^2\frac{m_1\frac{{\dot{X}}_{1,d}}{X_1}+m_2\frac{{\dot{X}}_{2,d}}{X_2}}{m_1+m_2}·\hfill \end{array}$$(E.4)

The temporal evolution of the variable 𝒳_j is hence simply given by 𝒳_j(t) = 𝒳_j(0)e^{(g_j + r_xj)t}. Orbits along the eccentric Lagrangian equilibria are defined by 𝒳₁ = 0, while those along the anti-Lagrangian equilibria are defined by 𝒳₂ = 0. For both of these configurations, the evolution of the eccentricities can be estimated from the quantity $X_j{\overline{X}}_j=e_j^2$:

$$\begin{array}{c}\hfill 2e_j{\dot{e}}_j={\dot{X}}_j{\overline{X}}_j+{\dot{\overline{X}}}_j{X}_j\end{array}$$(E.5)

where $\mathit{X}=P_{\mathcal{X}}^{-1}\mathcal{X}$, and $\dot{\mathit{X}}={\dot{P}}_{\mathcal{X}}^{-1}\mathcal{X}+P_{\mathcal{X}}^{-1}\dot{\mathcal{X}}$.

All Tables

All Figures