© ESO, 2017

1. Introduction

The angular momentum deficit (AMD)-stability criterion (Laskar 2000; Laskar & Petit 2017) allows to discriminate between a-priori stable planetary systems and systems needing an in-depth dynamical analysis to ensure their stability. The AMD-stability is based on the conservation of the angular momentum deficit (AMD, Laskar 1997) in the secular system at all orders of averaging (Laskar 2000; Laskar & Petit 2017). Indeed, the conservation of the AMD fixes an upper bound to the eccentricities. Since the semi-major axes are constant in the secular approximation, a low enough AMD forbids collisions between planets. The AMD-stability criterion has been used to classify planetary systems based on the stability of their secular dynamics (Laskar & Petit 2017).

However, while the analytical criterion developed in (Laskar & Petit 2017) does not depend on series expansions for small masses or spacing between the orbits, the secular hypothesis does not hold for systems experiencing mean motion resonances (MMR). Although a system with planets in MMR can be dynamically stable, chaotic behavior may result from the overlap of adjacent MMR, leading to a possible increase of the AMD and eventually to close encounters, collisions or ejections. For systems with small orbital separations, averaging over the mean anomalies is thus impossible due to the contribution of the first-order MMR terms. For example, two planets in circular orbits very close to each other are AMD-stable, however the dynamics of this system cannot be approximated by the secular dynamics. We thus need to modify the notion of AMD-stability in order to take into account those configurations.

In studies of planetary systems architecture, a minimal distance based on the Hill radius (Marchal & Bozis 1982) is often used as a criterion of stability (Gladman 1993; Chambers et al. 1996; Smith & Lissauer 2009; Pu & Wu 2015). However, Deck et al. (2013) suggested that stability criteria based on the MMR overlap are more accurate in characterizing the instability of the three-body planetary problem.

Based on the considerations of Chirikov (1979) for the overlap of resonant islands, Wisdom (1980) proposed a criterion of stability for the first-order MMR overlap in the context of the restricted circular three-body problem. This stability criterion defines a minimal distance between the orbits such that the first-order MMR overlap with one another. For orbits closer than this minimal distance, the MMR overlapping induces chaotic behavior eventually leading to the instability of the system.

Wisdom showed that the width of the chaotic region in the circular restricted problem is proportional to the ratio of the planet mass to the star mass to the power 2 / 7. Duncan et al. (1989) confirmed numerically that orbits closer than the Wisdom’s MMR overlap condition were indeed unstable. More recently, another stability criterion was proposed by Mustill & Wyatt (2012) to take into account the planet’s eccentricity. Deck et al. (2013) improved the two previous criteria by developing the resonant Hamiltonian for two massive, coplanar, low-eccentricity planets and Ramos et al. (2015) proposed a criterion of stability taking into account the second-order MMR in the restricted three-body problem.

While Deck’s criteria are in good agreement with numerical simulations (Deck et al. 2013) and can be applied to the three-body planetary problem, the case of circular orbits is still treated separately from the case of eccentric orbits. Indeed, the minimal distance imposed by the eccentric MMR overlap stability criterion vanishes with eccentricities and therefore cannot be applied to systems with small eccentricities. In this case, Mustill & Wyatt (2012) and Deck et al. (2013) use the criterion developed for circular orbits. A unified stability criterion for first-order MMR overlap had yet to be proposed.

In this paper, we propose in Sect. 2 a new derivation of the MMR overlap criterion based on the development of the three-body Hamiltonian by Delisle et al. (2012). We show in Sect. 3 how to obtain a unified criterion of stability working for both initially circular and eccentric orbits. In Sect. 4, we then use the defined stability criterion to limit the region where the dynamics can be considered to be secular and adapt the notion of AMD-stability thanks to the new limit of the secular dynamics. Finally we study in Sect. 5 how the modification of the AMD-stability definition affects the classification proposed in Laskar & Petit (2017).

2. The resonant Hamiltonian

The problem of two planets close to a first-order MMR on nearly circular and coplanar orbits can be reduced to a one-degree-of-freedom system through a sequence of canonical transformations (Wisdom 1986; Henrard et al. 1986; Delisle et al. 2012, 2014). We follow here the reduction of the Hamiltonian used in Delisle et al. (2012, 2014).

2.1. Averaged Hamiltonian in the vicinity of a resonance

Let us consider two planets of masses m₁ and m₂ orbiting a star of mass m₀ in the plane. We denote the positions of the planets, u_i, and the associated canonical momenta in the heliocentric frame, ${\begin{array}{c}{{˜}}_{}\\ {u}\end{array}}_{\mathit{i}}$. The Hamiltonian of the system is (Laskar & Robutel 1995) $\begin{array}{ccc}\mathrm{\mathscr{H}̂}& \mathrm{=}& \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathrm{2}}\left(\frac{\mathrm{1}}{\mathrm{2}}\frac{\mathrm{\parallel }{\begin{array}{c}{˜}\\ {u}\end{array}}_{\mathit{i}}{\mathrm{\parallel }}^{\mathrm{2}}}{{\mathit{m}}_{\mathit{i}}}\mathrm{-}\mathrm{𝒢}\frac{{\mathit{m}}_{\mathrm{0}}{\mathit{m}}_{\mathit{i}}}{{\mathit{u}}_{\mathit{i}}}\right)\\ & & \mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\frac{\mathrm{\parallel }{\begin{array}{c}{˜}\\ {u}\end{array}}_{\mathrm{1}}\mathrm{+}{\begin{array}{c}{˜}\\ {u}\end{array}}_{\mathrm{2}}{\mathrm{\parallel }}^{\mathrm{2}}}{{\mathit{m}}_{\mathrm{0}}}\mathrm{-}\mathrm{𝒢}\frac{{\mathit{m}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}}{{\mathrm{\Delta }}_{\mathrm{12}}}\end{array}$(1)where Δ₁₂ = ∥ u₁−u₂ ∥, and is the constant of gravitation. can be decomposed into a Keplerian part describing the motion of the planets if they had no masses and a perturbation part due to the influence of massive planets, $\begin{array}{ccc}\mathrm{𝒦̂}& \mathrm{=}& \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathrm{2}}\frac{\mathrm{1}}{\mathrm{2}}\frac{\mathrm{\parallel }{\begin{array}{c}{˜}\\ {u}\end{array}}_{\mathit{i}}{\mathrm{\parallel }}^{\mathrm{2}}}{{\mathit{m}}_{\mathit{i}}}\mathrm{-}\frac{\mathrm{𝒢}{\mathit{m}}_{\mathrm{0}}{\mathit{m}}_{\mathit{i}}}{{\mathit{u}}_{\mathit{i}}}\\ \mathit{\epsilon }\mathrm{\mathscr{H}̂}\mathrm{1}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{2}}\frac{\mathrm{\parallel }{\begin{array}{c}{˜}\\ {u}\end{array}}_{\mathrm{1}}\mathrm{+}{\begin{array}{c}{˜}\\ {u}\end{array}}_{\mathrm{2}}{\mathrm{\parallel }}^{\mathrm{2}}}{{\mathit{m}}_{\mathrm{0}}}\mathrm{-}\frac{\mathrm{𝒢}{\mathit{m}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}}{{\mathrm{\Delta }}_{\mathrm{12}}}\mathrm{·}\end{array}$The small parameter ε is defined as the ratio of the planet masses over the star mass $\mathit{\epsilon }\mathrm{=}\frac{{\mathit{m}}_{\mathrm{1}}\mathrm{+}{\mathit{m}}_{\mathrm{2}}}{{\mathit{m}}_{\mathrm{0}}}\mathrm{·}$(4)Let us denote the angular momentum, ${Ĝ}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathrm{2}}{{u}}_{\mathit{i}}\mathrm{\wedge }\begin{array}{c}{˜}\\ {u}\end{array}i$(5)which is simply the sum of the two planets Keplerian angular momentum. Ĝ is a first integral of the system.

Following Laskar (1991), we express the Hamiltonian in terms of the Poincaré coordinates $\begin{array}{ccc}\mathrm{\mathscr{H}̂}& \mathrm{=}& \mathrm{𝒦̂}\mathrm{+}\mathit{\epsilon }\mathrm{\mathscr{H}̂}\mathrm{1}\mathrm{\left(}\mathrm{\Lambda ̂}\mathit{i}\mathit{,}\mathit{x̂}\mathit{i}\mathit{,}\mathit{x̅}\mathit{i}\mathrm{\right)}\\ & \mathrm{=}& \mathrm{-}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathrm{2}}\frac{{\mathit{\mu }}^{\mathrm{2}}{\mathit{m}}_{\mathit{i}}^{\mathrm{3}}}{\mathrm{2}\mathrm{\Lambda ̂}\begin{array}{c}\mathrm{2}\\ \mathit{i}\end{array}}\mathrm{+}\mathit{\epsilon }\underset{\mathit{k}\mathrm{\in }{Z}^{\mathrm{2}}}{\sum _{\mathit{l,}\mathit{l̅}\mathrm{\in }{N}^{\mathrm{2}}}}{\mathit{C}}_{\mathit{l,}\mathit{l̅}\mathit{,k}}\mathrm{\left(}\mathrm{\Lambda ̂}\mathrm{\right)}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathrm{2}}{\mathrm{􏽙}}}\mathit{x̂}\begin{array}{c}{\mathit{l}}_{\mathit{i}}\\ \mathit{i}\end{array}\mathit{x̅}\begin{array}{c}\mathit{l̅}\mathit{i}\\ \mathit{i}\end{array}{\mathrm{e}}^{\mathrm{i}{\mathit{k}}_{\mathit{i}}{\mathit{\lambda }}_{\mathit{i}}}\mathit{,}\end{array}$(6)where and for i = 1,2, $\begin{array}{ccc}& & \mathrm{\Lambda ̂}\mathit{i}\mathrm{=}{\mathit{m}}_{\mathit{i}}\sqrt{\mathit{\mu }{\mathit{a}}_{\mathit{i}}}\\ & & \mathit{Ĝ}\mathit{i}\mathrm{=}\mathrm{\Lambda ̂}\mathit{i}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathit{i}}^{\mathrm{2}}}\\ & & \mathit{Ĉ}\mathit{i}\mathrm{=}\mathrm{\Lambda ̂}\mathit{i}\mathrm{-}\mathit{Ĝ}\mathit{i}\\ & & \mathit{x̂}\mathit{i}\mathrm{=}\sqrt{\mathit{Ĉ}\mathit{i}}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{\mathit{\varpi }}_{\mathit{i}}}\\ & & {\mathit{\lambda }}_{\mathit{i}}\mathrm{=}{\mathit{M}}_{\mathit{i}}\mathrm{+}{\mathit{\varpi }}_{\mathit{i}}\mathit{.}\end{array}$Here, M_i corresponds to the mean anomaly, ϖ_i to the longitude of the pericenter, a_i to the semi-major axis and e_i to the eccentricity of the Keplerian orbit of the planet i. Ĝ_i is the Keplerian angular momentum of planet i. We use the set of symplectic coordinates of the problem , or the canonically associated variables . The coefficients depend on and the masses of the bodies. They are linear combinations of Laplace coefficients (Laskar & Robutel 1995). As a consequence of angular momentum conservation, the d’Alembert rule gives a relation on the indices of the non-zero coefficients $\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathrm{2}}{\mathit{k}}_{\mathit{i}}\mathrm{-}{\mathit{l}}_{\mathit{i}}\mathrm{+}\mathit{l̅}\mathit{i}\mathrm{=}\mathrm{0.}$(7)We study here a system with periods close to the first-order MMR p:p + 1 with p ∈ N^∗. For periods close to this configuration, we have −pn₁ + (p + 1)n₂ ≃ 0, where ${\mathit{n}}_{\mathit{i}}\mathrm{=}{\mathit{\mu }}^{\mathrm{2}}{\mathit{m}}_{\mathit{i}}^{\mathrm{3}}\mathit{/}\mathrm{\Lambda ̂}\begin{array}{c}\mathrm{3}\\ \mathit{i}\end{array}$ is the Keplerian mean motion of the planet i.

2.1.1. Averaging over non-resonant mean-motions

Due to the p:p + 1 resonance, we cannot average on both mean anomalies independently. Therefore, there is no conservation of as in the secular problem. However, the partial averaging over one of the mean anomaly gives another first integral. Following (Delisle et al. 2012), we consider the equivalent set of coordinates , and make the following change of angles $\left(\begin{array}{c}\\ \mathit{\sigma }\\ {\mathit{M}}_{\mathrm{2}}\end{array}\right)\mathrm{=}\left(\begin{array}{c}\\ \mathrm{-}\mathit{p}& \mathit{p}\mathrm{+}\mathrm{1}\\ \mathrm{0}& \mathrm{1}\end{array}\right)\left(\begin{array}{c}\\ {\mathit{M}}_{\mathrm{1}}\\ {\mathit{M}}_{\mathrm{2}}\end{array}\right)\mathit{.}$(8)The actions associated to these angles are $\left(\begin{array}{c}\\ \mathrm{\Gamma ̂}\mathrm{1}\\ \mathrm{\Gamma ̂}\end{array}\right)\mathrm{=}\left(\begin{array}{c}\\ \mathrm{-}\frac{\mathrm{1}}{\mathit{p}}& \mathrm{0}\\ \frac{\mathit{p}\mathrm{+}\mathrm{1}}{\mathit{p}}& \mathrm{1}\end{array}\right)\left(\begin{array}{c}\\ \mathrm{\Lambda ̂}\mathrm{1}\\ \mathrm{\Lambda ̂}\mathrm{2}\end{array}\right)\mathrm{=}\left(\begin{array}{c}\\ \frac{-1}{\mathit{p}}\mathrm{\Lambda ̂}\mathrm{1}\\ \frac{\mathit{p}\mathrm{+}\mathrm{1}}{\mathit{p}}\mathrm{\Lambda ̂}\mathrm{1}\mathrm{+}\mathrm{\Lambda ̂}\mathrm{2}\end{array}\right)\mathit{.}$(9)We can now average the Hamiltonian over M₂ using a change of variables close to the identity given by the Lie series method. Up to terms of orders ε², we can kill all the terms with indices not of the form . In order to keep the notations light, we do not change the name of the variables after the averaging. We also designate the remaining coefficients by the lighter expression . Since M₂ does not appear explicitly in the remaining terms, $\mathrm{\Gamma ̂}\mathrm{=}\frac{\mathit{p}\mathrm{+}\mathrm{1}}{\mathit{p}}\mathrm{\Lambda ̂}\mathrm{1}\mathrm{+}\mathrm{\Lambda ̂}\mathrm{2}$(10)is a first integral of the averaged Hamiltonian. The parameter is often designed as the spacing parameter (Michtchenko et al. 2008) and has been used extensively in the study of the first-order MMR dynamics.

Expressed with the variables , the Hamiltonian can be written $\begin{array}{ccc}\mathrm{\mathscr{H}̂}\mathrm{av}& \mathrm{=}& \mathrm{-}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathrm{2}}\frac{{\mathit{\mu }}^{\mathrm{2}}{\mathit{m}}_{\mathit{i}}^{\mathrm{3}}}{\mathrm{2}\mathrm{\Lambda ̂}\begin{array}{c}\mathrm{2}\\ \mathit{i}\end{array}}\\ & & \end{array}$(11)where we dropped the terms of order ε².

2.1.2. Poincare-like complex coordinates

Delisle et al. (2012) used a change of the angular coordinates in order to remove the exponential in the second term of Eq. (11) and use Ĝ and as actions. The new set of angles (θ_Γ,θ_G,σ₁,σ₂) is defined as $\left(\begin{array}{c}\\ {\mathit{\theta }}_{\mathrm{\Gamma }}\\ {\mathit{\theta }}_{\mathit{G}}\\ {\mathit{\sigma }}_{\mathrm{1}}\\ {\mathit{\sigma }}_{\mathrm{2}}\end{array}\right)\mathrm{=}\left(\begin{array}{c}\\ \mathit{p}& \mathrm{-}\mathit{p}& \mathrm{0}& \mathrm{0}\\ \mathrm{-}\mathit{p}& \mathit{p}\mathrm{+}\mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{-}\mathit{p}& \mathit{p}\mathrm{+}\mathrm{1}& \mathrm{1}& \mathrm{0}\\ \mathrm{-}\mathit{p}& \mathit{p}\mathrm{+}\mathrm{1}& \mathrm{0}& \mathrm{1}\end{array}\right)\mathrm{·}\left(\begin{array}{c}\\ {\mathit{\lambda }}_{\mathrm{1}}\\ {\mathit{\lambda }}_{\mathrm{2}}\\ \mathrm{-}{\mathit{\varpi }}_{\mathrm{1}}\\ \mathrm{-}{\mathit{\varpi }}_{\mathrm{2}}\end{array}\right)\mathit{.}$(12)The conjugated actions are $\left(\begin{array}{c}\\ \mathrm{\Gamma ̂}\\ \mathit{Ĝ}\\ \mathit{Ĉ}\mathrm{1}\\ \mathit{Ĉ}\mathrm{2}\end{array}\right)\mathrm{=}\left(\begin{array}{c}\\ \frac{\mathit{p}\mathrm{+}\mathrm{1}}{\mathit{p}}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}& \mathrm{1}& -1& -1\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right)\mathrm{·}\left(\begin{array}{c}\\ \mathrm{\Lambda ̂}\mathrm{1}\\ \mathrm{\Lambda ̂}\mathrm{2}\\ \mathit{Ĉ}\mathrm{1}\\ \mathit{Ĉ}\mathrm{2}\end{array}\right)\mathit{.}$(13)We define $\hat{X}{}_{\mathit{i}}\mathrm{=}\sqrt{\mathit{Ĉ}{}_{\mathit{i}}}{\mathrm{e}}^{\mathrm{i}{\mathit{\sigma }}_{\mathit{i}}}$, the complex coordinates associated to (Ĉ_i,σ_i). Since we have $\hat{X}{}_{\mathit{i}}\mathrm{=}\mathit{x̂}{}_{\mathit{i}}\mathrm{e}^{\mathrm{i}{\mathit{\theta }}_{\mathit{G}}}$, the terms of the perturbation in Eq. (11) can be written $\begin{array}{ccc}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathrm{2}}{\mathrm{􏽙}}}\mathit{x̂}\begin{array}{c}{\mathit{l}}_{\mathit{i}}\\ \mathit{i}\end{array}\mathit{x̅}\begin{array}{c}\mathit{l̅}\mathit{i}\\ \mathit{i}\end{array}{\mathrm{e}}^{\mathrm{i}\mathit{j}{\mathit{\theta }}_{\mathit{G}}}& \mathrm{=}& \underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathrm{2}}{\mathrm{􏽙}}}\hat{X}\begin{array}{c}{\mathit{l}}_{\mathit{i}}\\ \mathit{i}\end{array}\overline{X}\begin{array}{c}\mathit{l̅}\mathit{i}\\ \mathit{i}\end{array}{\mathrm{e}}^{\mathrm{i}\mathrm{(}\mathrm{-}{\mathit{l}}_{\mathit{i}}\mathrm{+}\mathit{l̅}\mathit{i}\mathrm{+}\mathit{j}\mathrm{)}{\mathit{\theta }}_{\mathit{G}}}\\ & \mathrm{=}& \underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathrm{2}}{\mathrm{􏽙}}}\hat{X}\begin{array}{c}{\mathit{l}}_{\mathit{i}}\\ \mathit{i}\end{array}\overline{X}\begin{array}{c}\mathit{l̅}\mathit{i}\\ \mathit{i}\end{array}\mathrm{;}\end{array}$(14)the last equality resulting from the d’Alembert rule (7). and Ĝ are conserved and the averaged Hamiltonian no longer depends on the angles θ_Γ and θ_G$\mathrm{\mathscr{H}̂}\mathrm{av}\mathrm{=}\mathrm{-}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathrm{2}}\frac{{\mathit{\mu }}^{\mathrm{2}}{\mathit{m}}_{\mathit{i}}^{\mathrm{3}}}{\mathrm{2}\mathrm{\Lambda ̂}\begin{array}{c}\mathrm{2}\\ \mathit{i}\end{array}}\mathrm{+}\mathit{\epsilon }\underset{\mathit{j}\mathrm{\in }Z}{\sum _{\mathit{l,}\mathit{l̅}\mathrm{\in }{N}^{\mathrm{2}}}}{\mathit{C}}_{\mathit{l,}\mathit{l̅}\mathit{,j}}\mathrm{\left(}\mathrm{\Lambda ̂}\mathrm{\right)}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathrm{2}}{\mathrm{􏽙}}}\hat{X}\begin{array}{c}{\mathit{l}}_{\mathit{i}}\\ \mathit{i}\end{array}\overline{X}\begin{array}{c}\mathit{l̅}\mathit{i}\\ \mathit{i}\end{array}\mathit{.}$(15)and can be expressed as functions of the new variables and we have $\begin{array}{ccc}\mathrm{\Lambda ̂}\mathrm{1}& \mathrm{=}& \mathrm{-}\mathit{p}\mathrm{(}\mathit{Ĉ}\mathrm{+}\mathit{Ĝ}\mathrm{-}\mathrm{\Gamma ̂}\mathrm{)}\\ \mathrm{\Lambda ̂}\mathrm{2}& \mathrm{=}& \mathrm{(}\mathit{p}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{(}\mathit{Ĉ}\mathrm{+}\mathit{Ĝ}\mathrm{)}\mathrm{-}\mathit{p}\mathrm{\Gamma ̂}\mathit{,}\end{array}$where Ĉ = Ĉ₁ + Ĉ₂ is the total AMD of the system. Up to the value of the first integrals and Ĝ, the system now has two effective degrees of freedom.

2.2. Computation of the perturbation coefficients

We now truncate the perturbation, keeping only the leading-order terms. Since we consider the first-order MMR, the Hamiltonian contains some linear terms in X_i. Therefore the secular terms are neglected since they are at least quadratic. Moreover, the restriction to the planar problem is justified since the inclination terms are at least of order two.

We follow the method described in Laskar (1991) and Laskar & Robutel (1995) to determine the expression of the perturbation . The details of the computation are given in Appendix A. Since we compute an expression at first order in eccentricities and ε, the semi major axis and in particular their ratio, $\mathit{\alpha }\mathrm{=}\frac{{\mathit{a}}_{\mathrm{1}}}{{\mathit{a}}_{\mathrm{2}}}\mathit{,}$(18)are evaluated at the resonance. At the first order, the perturbation term has for expression $\mathit{\epsilon }\mathrm{\mathscr{H}̂}\mathrm{1}\mathrm{=}\mathit{R̂}\mathrm{1}\mathrm{(}\hat{X}\mathrm{1}\mathrm{+}\overline{X}\mathrm{1}\mathrm{)}\mathrm{+}\mathit{R̂}\mathrm{2}\mathrm{(}\hat{X}\mathrm{2}\mathrm{+}\overline{X}\mathrm{2}\mathrm{)}\mathit{,}$(19)where $\mathit{R̂}\mathrm{1}\mathrm{=}\mathrm{-}\mathit{\epsilon }\frac{\mathit{\gamma }}{\mathrm{1}\mathrm{+}\mathit{\gamma }}\frac{{\mathit{\mu }}^{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}^{\mathrm{3}}}{\mathrm{\Lambda ̂}\begin{array}{c}\mathrm{2}\\ \mathrm{2}\end{array}}\frac{\mathrm{1}}{\mathrm{2}}\sqrt{\frac{\mathrm{2}}{\mathrm{\Lambda ̂}\mathrm{1}}}{\mathit{r}}_{\mathrm{1}}\mathrm{\left(}\mathit{\alpha }\mathrm{\right)}$(20)and $\mathit{R̂}\mathrm{2}\mathrm{=}\mathrm{-}\mathit{\epsilon }\frac{\mathit{\gamma }}{\mathrm{1}\mathrm{+}\mathit{\gamma }}\frac{{\mathit{\mu }}^{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}^{\mathrm{3}}}{\mathrm{\Lambda ̂}\begin{array}{c}\mathrm{2}\\ \mathrm{2}\end{array}}\frac{\mathrm{1}}{\mathrm{2}}\sqrt{\frac{\mathrm{2}}{\mathrm{\Lambda ̂}\mathrm{2}}}{\mathit{r}}_{\mathrm{2}}\mathrm{\left(}\mathit{\alpha }\mathrm{\right)}$(21)with γ = m₁/m₂, $\begin{array}{ccc}{\mathit{r}}_{\mathrm{1}}\mathrm{\left(}\mathit{\alpha }\mathrm{\right)}& \mathrm{=}& \\ \mathrm{and}& & \\ {\mathit{r}}_{\mathrm{2}}\mathrm{\left(}\mathit{\alpha }\mathrm{\right)}& \mathrm{=}& \frac{\mathit{\alpha }}{\mathrm{4}}\mathrm{(}\mathrm{3}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{(}\mathit{p}\mathrm{-}\mathrm{1}\mathrm{)}}\mathrm{\left(}\mathit{\alpha }\mathrm{\right)}\mathrm{-}\mathrm{2}\mathit{\alpha }{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathit{p}\mathrm{\right)}}\mathrm{\left(}\mathit{\alpha }\mathrm{\right)}\mathrm{-}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{(}\mathit{p}\mathrm{+}\mathrm{1}\mathrm{)}}\mathrm{(}\mathit{\alpha }{\mathrm{)}}^{\mathrm{)}}\\ & & \end{array}$In the two previous expressions, ${\mathit{b}}_{\mathit{s}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}\mathrm{\left(}\mathit{\alpha }\mathrm{\right)}$ are the Laplace coefficients that can be expressed as ${\mathit{b}}_{\mathit{s}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}\mathrm{\left(}\mathit{\alpha }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathit{\pi }}{\mathrm{\int }}_{\mathrm{-}\mathit{\pi }}^{\mathit{\pi }}\frac{\cos \mathrm{\left(}\mathit{k\phi }\mathrm{\right)}}{{\left(\mathrm{1}\mathrm{-}\mathrm{2}\mathit{\alpha }\cos \mathit{\phi }\mathrm{+}{\mathit{\alpha }}^{\mathrm{2}}\right)}^{\mathit{s}}}\mathrm{d}\mathit{\phi }$(24)for k> 0. For k = 0, a 1 / 2 factor has to be added in the second-hand member of (24).

For p = 1, it should be noted that a contribution from the kinetic part should be added (Appendix A and Delisle et al. 2012) ${\mathrm{\mathscr{H}}}_{\mathrm{1}\mathit{,i}}\mathrm{=}\frac{{\mathit{\mu }}^{\mathrm{2}}{\mathit{m}}_{\mathrm{1}}^{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}^{\mathrm{2}}}{\mathrm{2}{\mathit{m}}_{\mathrm{0}}\mathrm{\Lambda ̂}\mathrm{1}\mathrm{\Lambda ̂}\mathrm{2}}\sqrt{\frac{\mathrm{2}}{\mathrm{\Lambda ̂}\mathrm{2}}}\mathrm{(}\hat{X}\mathrm{2}\mathrm{+}\overline{X}\mathrm{2}\mathrm{)}\mathit{.}$(25)Using the expression of α at the resonance p:p + 1, ${\mathit{\alpha }}_{\mathrm{0}}\mathrm{=}{\left(\frac{\mathit{p}}{\mathit{p}\mathrm{+}\mathrm{1}}\right)}^{\mathrm{2}\mathit{/}\mathrm{3}}\mathit{,}$(26)we can give the asymptotic development of the coefficients r₁ and r₂ for p → + ∞ (see Appendix A.1). The equivalent is $\mathrm{-}{\mathit{r}}_{\mathrm{1}}\mathrm{~}{\mathit{r}}_{\mathrm{2}}\mathrm{~}\frac{{\mathit{K}}_{\mathrm{1}}\mathrm{(}\mathrm{2}\mathit{/}\mathrm{3}\mathrm{)}\mathrm{+}\mathrm{2}{\mathit{K}}_{\mathrm{0}}\mathrm{(}\mathrm{2}\mathit{/}\mathrm{3}\mathrm{)}}{\mathit{\pi }}\mathrm{(}\mathit{p}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{;}$(27)where K_ν(x) is the modified Bessel function of the second kind. We note r the numerical factor of the equivalent (27), we have $\mathit{r}\mathrm{=}\frac{{\mathit{K}}_{\mathrm{1}}\mathrm{(}\mathrm{2}\mathit{/}\mathrm{3}\mathrm{)}\mathrm{+}\mathrm{2}{\mathit{K}}_{\mathrm{0}}\mathrm{(}\mathrm{2}\mathit{/}\mathrm{3}\mathrm{)}}{\mathit{\pi }}\mathrm{=}\mathrm{0.80199.}$(28)For the resonant coefficients r₁ and r₂, Deck et al. (2013) used the expressions f_{p + 1,27}(α) and f_{p + 1,31}(α) given in Murray & Dermott (1999, pp. 539−556). The expressions (22) and (23) are similar to f_{p + 1,27}(α) and f_{p + 1,31}(α) up to algebraic transformations using the relations between Laplace coefficients (Laskar & Robutel 1995). In their computations, Deck et al. used a numerical fit of the coefficients for p = 2 to 150 and obtained $\mathrm{-}{\mathit{f}}_{\mathit{p}\mathrm{+}\mathrm{1}\mathit{,}\mathrm{27}}\mathrm{~}{\mathit{f}}_{\mathit{p}\mathrm{+}\mathrm{1}\mathit{,}\mathrm{31}}\mathrm{~}\mathrm{0.802}\mathit{p}\mathit{.}$(29)We obtain the same numerical factor r through the analytical development of the functions r₁ and r₂.

2.3. Renormalization

So far, the Hamiltonian has two degrees of freedom $\mathrm{\left(}\hat{X}{}_{\mathrm{1}}\mathit{,}\overline{X}{}_{\mathrm{1}}\mathit{,}\hat{X}{}_{\mathrm{2}}\mathit{,}\overline{X}{}_{\mathrm{2}}\mathrm{\right)}$ and depends on two parameters Ĝ and . As shown in Delisle et al. (2012), the constant can be used to scale the actions, the Hamiltonian and the time without modifying the dynamics. We define $\begin{array}{ccc}& & {\mathrm{\Lambda }}_{\mathit{i}}\mathrm{=}\mathrm{\Lambda ̂}\mathit{i}\mathit{/}\mathrm{\Gamma ̂}\mathit{,}\\ & & \mathit{G}\mathrm{=}\mathit{Ĝ}\mathit{/}\mathrm{\Gamma ̂}\mathit{,}\\ & & {\mathit{C}}_{\mathit{i}}\mathrm{=}\mathit{Ĉ}\mathit{i}\mathit{/}\mathrm{\Gamma ̂}\mathit{,}\\ & & X\mathit{i}\mathrm{=}\hat{X}\mathit{i}\mathit{/}\sqrt{\mathrm{\Gamma ̂}}\mathit{,}\\ & & \mathrm{\mathscr{H}}\mathrm{=}\mathrm{\Gamma ̂}\mathrm{2}\mathrm{\mathscr{H}̂}\mathit{,}\\ & & \mathit{t}\mathrm{=}\mathit{t̂}\mathit{/}\mathrm{\Gamma ̂}\mathrm{3}\mathit{.}\end{array}$With this change of variables, the new Hamiltonian no longer depends on .

The shape of the phase space is now only dependent on the first integral G. However, G does not vanish for the configuration around which the Hamiltonian is developed: the case of two resonant planets on circular orbits. To be able to develop the Keplerian part in power of the system’s parameter, we define ΔG = G₀−G, the difference in angular momentum between the circular resonant system and the actual configuration. We have ${\mathit{G}}_{\mathrm{0}}\mathrm{=}{\mathrm{\Lambda }}_{\mathrm{1}\mathit{,}\mathrm{0}}\mathrm{+}{\mathrm{\Lambda }}_{\mathrm{2}\mathit{,}\mathrm{0}}\mathit{,}$(30)where Λ_1,0 and Λ_2,0 are the value of Λ₁ and Λ₂ at resonance. By definition, we have $\frac{{\mathrm{\Lambda }}_{\mathrm{1}\mathit{,}\mathrm{0}}}{{\mathrm{\Lambda }}_{\mathrm{2}\mathit{,}\mathrm{0}}}\mathrm{=}\mathit{\gamma }{\left(\frac{\mathit{p}}{\mathit{p}\mathrm{+}\mathrm{1}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathrm{=}\mathit{\gamma }\sqrt{{\mathit{\alpha }}_{\mathrm{0}}}\mathit{.}$(31)Moreover, we can express Λ_1,0 as a function of the ratios α₀ and γ, ${\mathrm{\Lambda }}_{\mathrm{1}\mathit{,}\mathrm{0}}\mathrm{=}\frac{\mathrm{\Lambda ̂}\mathrm{1}\mathit{,}\mathrm{0}}{\mathrm{\Gamma ̂}\mathrm{0}}\mathrm{=}\frac{\mathrm{1}}{\left(\frac{\mathit{p}\mathrm{+}\mathrm{1}}{\mathit{p}}\right)\mathrm{+}\frac{{\mathrm{\Lambda }}_{\mathrm{2}\mathit{,}\mathrm{0}}}{{\mathrm{\Lambda }}_{\mathrm{1}\mathit{,}\mathrm{0}}}}\mathrm{=}\left(\frac{\mathit{p}}{\mathit{p}\mathrm{+}\mathrm{1}}\right)\frac{\mathit{\gamma }}{\mathit{\gamma }\mathrm{+}{\mathit{\alpha }}_{\mathrm{0}}}\mathrm{·}$(32)Similarly, Λ_2,0 can be expressed as ${\mathrm{\Lambda }}_{\mathrm{2}\mathit{,}\mathrm{0}}\mathrm{=}\frac{{\mathit{\alpha }}_{\mathrm{0}}}{{\mathit{\alpha }}_{\mathrm{0}}\mathrm{+}\mathit{\gamma }}\mathrm{·}$(33)Since G₀ is constant, ΔG is also a first integral of ℋ. From now on, we consider ΔG as a parameter of the two-degrees-of-freedom (X₁,X₂) Hamiltonian ℋ. The Keplerian part depends on the coordinates X_i through the dependence of Λ_i in C.

Λ₁and Λ₂ can be expressed as functions of the Hamiltonian coordinates and their value at the resonance, $\begin{array}{ccc}{\mathrm{\Lambda }}_{\mathrm{1}}& \mathrm{=}& {\mathrm{\Lambda }}_{\mathrm{1}\mathit{,}\mathrm{0}}\mathrm{-}\mathit{p}\mathrm{(}\mathit{C}\mathrm{-}\mathrm{\Delta }\mathit{G}\mathrm{)}\\ {\mathrm{\Lambda }}_{\mathrm{2}}& \mathrm{=}& \end{array}$(34)

2.4. Integrable Hamiltonian

The system can be made integrable by a rotation of the coordinates X_i (Sessin & Ferraz-Mello 1984; Henrard et al. 1986; Delisle et al. 2014). We introduce R and φ such that ${\mathit{R}}_{\mathrm{1}}\mathrm{=}\mathit{R}\cos \mathrm{\left(}\mathit{\phi }\mathrm{\right)} \mathrm{and} {\mathit{R}}_{\mathrm{2}}\mathrm{=}\mathit{R}\sin \mathrm{\left(}\mathit{\phi }\mathrm{\right)}\mathit{.}$(35)We have ${\mathit{R}}^{\mathrm{2}}\mathrm{=}{\mathit{R}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{\mathit{R}}_{\mathrm{2}}^{\mathrm{2}}$ and tan(φ) = R₂/R₁. If we note ℛ_φ the rotation of angle φ we define y such that X= ℛ_φy. We still have so the only change in the Hamiltonian is the perturbation term $\begin{array}{ccc}\mathrm{\mathscr{H}}& \mathrm{=}& \mathrm{𝒦}\mathrm{\left(}\mathit{C,}\mathrm{\Delta }\mathit{G}\mathrm{\right)}\mathrm{+}\mathit{R}\mathrm{(}{\mathit{y}}_{\mathrm{1}}\mathrm{+}\mathit{y̅}\mathrm{1}\mathrm{)}\\ & \mathrm{=}& \mathrm{𝒦}\mathrm{\left(}\mathit{C,}\mathrm{\Delta }\mathit{G}\mathrm{\right)}\mathrm{+}\mathrm{2}\mathit{R}\sqrt{{\mathit{I}}_{\mathrm{1}}}\cos \mathrm{\left(}{\mathit{\theta }}_{\mathrm{1}}\mathrm{\right)}\mathit{,}\end{array}$(36)where (I,θ) are the action-angle coordinates associated to y. With these coordinates, I₂ is a first integral. R has for expression ${\mathit{R}}^{\mathrm{2}}\mathrm{=}{\left(\frac{\mathit{\epsilon \gamma }}{\mathrm{1}\mathrm{+}\mathit{\gamma }}\frac{{\mathit{\mu }}^{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}^{\mathrm{3}}}{{\mathrm{\Lambda }}_{\mathrm{2}\mathit{,}\mathrm{0}}^{\mathrm{2}}}\right)}^{\mathrm{2}}\left(\frac{{\mathit{r}}_{\mathrm{1}}\mathrm{(}{\mathit{\alpha }}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{2}{\mathrm{\Lambda }}_{\mathrm{1}\mathit{,}\mathrm{0}}}\mathrm{+}\frac{{\mathit{r}}_{\mathrm{2}}\mathrm{(}{\mathit{\alpha }}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{2}{\mathrm{\Lambda }}_{\mathrm{2}\mathit{,}\mathrm{0}}}\right)\mathrm{·}$(37)We now develop the Keplerian part around the circular resonant configuration in series of (C−ΔG) thanks to the relations (34). We develop the Keplerian part to the second order in (C−ΔG) since the first order vanishes (see Appendix B). The computation of the second-order coefficient gives $\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{𝒦}}_{\mathrm{2}}\mathrm{=}\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}{\mathit{\mu }}^{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}^{\mathrm{3}}\frac{\mathrm{(}\mathit{\gamma }\mathrm{+}{\mathit{\alpha }}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{5}}}{\mathit{\gamma }{\mathit{\alpha }}_{\mathrm{0}}^{\mathrm{4}}}\mathrm{(}\mathit{p}\mathrm{+}\mathrm{1}{\mathrm{)}}^{\mathrm{2}}\mathit{.}$(38)We drop the constant part of the Hamiltonian and obtain the following expression $\mathrm{\mathscr{H}}\mathrm{=}\frac{{\mathrm{𝒦}}_{\mathrm{2}}}{\mathrm{2}}\mathrm{(}{\mathit{I}}_{\mathrm{1}}\mathrm{+}{\mathit{I}}_{\mathrm{2}}\mathrm{-}\mathrm{\Delta }\mathit{G}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{2}\mathit{R}\sqrt{{\mathit{I}}_{\mathrm{1}}}\cos \mathrm{(}{\mathit{\theta }}_{\mathrm{1}}\mathrm{)}\mathit{.}$(39)We again change the time scale by dividing the Hamiltonian by and multiplying the time by this factor. We define $\mathit{\chi }\mathrm{=}\mathrm{-}\frac{\sqrt{\mathrm{2}}\mathit{R}}{{\mathrm{𝒦}}_{\mathrm{2}}}$(40)and after simplification, $\begin{array}{ccc}\mathit{\chi }& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{3}}\frac{\mathit{\epsilon }\mathrm{(}\mathit{\gamma }{\mathit{\alpha }}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\gamma }\mathrm{)}\mathrm{(}{\mathit{\alpha }}_{\mathrm{0}}\mathrm{+}\mathit{\gamma }{\mathrm{)}}^{\mathrm{2}}}\frac{{\mathit{r}}_{\mathrm{2}}\mathrm{\left(}{\mathit{\alpha }}_{\mathrm{0}}\mathrm{\right)}}{\mathrm{(}\mathit{p}\mathrm{+}\mathrm{1}{\mathrm{)}}^{\mathrm{2}}}\mathit{f}\mathrm{\left(}\mathit{p}\mathrm{\right)}\\ & \mathrm{=}& \frac{\mathit{r}}{\mathrm{3}}\frac{\mathit{\epsilon }{\mathit{\gamma }}^{\mathrm{3}\mathit{/}\mathrm{2}}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\gamma }{\mathrm{)}}^{\mathrm{3}}}\frac{\mathrm{1}}{\mathit{p}\mathrm{+}\mathrm{1}}\mathrm{+}\mathrm{O}\mathrm{\left(}\mathrm{\right(}\mathit{p}\mathrm{+}\mathrm{1}{\mathrm{)}}^{-2}\mathrm{)}\mathit{,}\end{array}$where r was defined in (28) and f(p) = 1 + O(p^-1) is a function of p and γ$\mathit{f}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{=}\sqrt{\mathrm{1}\mathrm{-}\frac{{\mathit{\alpha }}_{\mathrm{0}}}{{\mathit{\alpha }}_{\mathrm{0}}\mathrm{+}\mathit{\gamma }}\left(\mathrm{1}\mathrm{-}\frac{\mathit{p}\mathrm{+}\mathrm{1}}{\mathit{p}}{\left(\frac{{\mathit{r}}_{\mathrm{1}}}{{\mathit{r}}_{\mathrm{2}}}\right)}^{\mathrm{2}}\right)}\mathrm{·}$(43)At this point the Hamiltonian can be written $\mathrm{\mathscr{H}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}{\mathit{I}}_{\mathrm{1}}\mathrm{+}{\mathit{I}}_{\mathrm{2}}\mathrm{-}\mathrm{\Delta }\mathit{G}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathit{\chi }\sqrt{\mathrm{2}{\mathit{I}}_{\mathrm{1}}}\cos \mathrm{(}{\mathit{\theta }}_{\mathrm{1}}\mathrm{)}$(44)and has almost its final form. We divide the actions and the time by χ^{2 / 3} and the Hamiltonian by χ^{4 / 3} and we obtain ${\mathrm{\mathscr{H}}}_{\mathrm{A}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}\mathrm{ℐ}\mathrm{-}{\mathrm{ℐ}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}\mathrm{-}\sqrt{\mathrm{2}\mathrm{ℐ}}\cos \mathrm{(}{\mathit{\theta }}_{\mathrm{1}}\mathrm{)}\mathit{,}$(45)where $\mathrm{ℐ}\mathrm{=}{\mathit{\chi }}^{\mathrm{-}\mathrm{2}\mathit{/}\mathrm{3}}{\mathit{I}}_{\mathrm{1}} \mathrm{and} {\mathrm{ℐ}}_{\mathrm{0}}\mathrm{=}{\mathit{\chi }}^{\mathrm{-}\mathrm{2}\mathit{/}\mathrm{3}}\mathrm{(}\mathrm{\Delta }\mathit{G}\mathrm{-}{\mathit{I}}_{\mathrm{2}}\mathrm{)}\mathit{.}$(46)

2.5. Andoyer Hamiltonian

We now perform a polar to Cartesian change of coordinates with $\begin{array}{ccc}\mathit{X}& \mathrm{=}& \mathrm{-}\sqrt{\mathrm{2}\mathrm{ℐ}}\cos \mathrm{\left(}{\mathit{\theta }}_{\mathrm{1}}\mathrm{\right)}\mathit{,}\\ \mathit{Y}& \mathrm{=}& \end{array}$(47)We change the sign of X in order to have the same orientation as (Deck et al. 2013). Doing so, the Hamiltonian becomes ${\mathrm{\mathscr{H}}}_{\mathrm{A}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}{\mathit{X}}^{\mathrm{2}}\mathrm{+}{\mathit{Y}}^{\mathrm{2}}\mathrm{)}\mathrm{-}{\mathrm{ℐ}}_{\mathrm{0}}\right)}^{\mathrm{2}}\mathrm{-}\mathit{X}\mathit{.}$(48)We recognize the second fundamental model of resonance (Henrard & Lemaitre 1983). This Hamiltonian is also called an Andoyer Hamiltonian (Ferraz-Mello 2007). We show in Fig. 1 the level curves of the Hamiltonian ℋ_A for ℐ₀ = 3.

The fixed points of the Hamiltonian satisfy the equations $\begin{array}{ccc}\mathit{Ẋ}& \mathrm{=}& \mathit{Y}\left(\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}{\mathit{X}}^{\mathrm{2}}\mathrm{+}{\mathit{Y}}^{\mathrm{2}}\mathrm{)}\mathrm{-}{\mathrm{ℐ}}_{\mathrm{0}}\right)\mathrm{=}\mathrm{0}\\ \mathit{Ẏ}& \mathrm{=}& \mathrm{-}\mathit{X}\left(\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}{\mathit{X}}^{\mathrm{2}}\mathrm{+}{\mathit{Y}}^{\mathrm{2}}\mathrm{)}\mathrm{-}{\mathrm{ℐ}}_{\mathrm{0}}\right)\mathrm{-}\mathrm{1}\mathrm{=}\mathrm{0}\mathit{,}\end{array}$which have for solutions Y = 0 and the real roots of the cubic equation in X${\mathit{X}}^{\mathrm{3}}\mathrm{-}\mathrm{2}{\mathrm{ℐ}}_{\mathrm{0}}\mathit{X}\mathrm{+}\mathrm{2}\mathrm{=}\mathrm{0.}$(51)Equation (51) has three solutions (Deck et al. 2013) if its determinant $\mathrm{\Delta }\mathrm{=}\mathrm{32}\mathrm{(}{\mathrm{ℐ}}_{\mathrm{0}}^{\mathrm{3}}\mathrm{-}\mathrm{27}\mathit{/}\mathrm{8}\mathrm{)}\mathit{>}\mathrm{0}$, i.e. ℐ₀> 3 / 2. In this case, we note these roots X₁<X₂<X₃. X₁ and X₂ are elliptic fixed points while X₃ is a hyperbolic one.

Fig. 1 Hamiltonian ℋA (48) represented with the saddle point and the separatrices in red.

3. Overlap criterion

As seen in the previous section, the motion of two planets near a first-order MMR can be reduced to an integrable system for small eccentricities and planet masses. However, if two independent combinations of frequencies are close to zero at the same time, the previous reduction is not valid anymore. Indeed, we must then keep, in the averaging, the terms corresponding to both resonances. While for a single resonant term the system is integrable, overlapping resonant islands will lead to chaotic motion (Chirikov 1979).

Wisdom (1980) first applied the resonance overlap criterion to the first-order MMR and found, in the case of the restricted three-body problem with a circular planet, that the overlap occurs for $\mathrm{1}\mathrm{-}\mathit{\alpha }\mathit{<}\mathrm{1.3}{\mathit{\epsilon }}^{\mathrm{2}\mathit{/}\mathrm{7}}\mathit{.}$(52)Through numerical simulations, (Duncan et al. 1989) confirmed Wisdom’s expression up to the numerical coefficient {(1−α< 1.5 ε^{2 / 7})}. A similar criterion was then developed by Mustill & Wyatt (2012) for an eccentric planet, they found that for an eccentricity above 0.2 ε^{3 / 7}, the overlap region satisfies the criterion 1−α< 1.8(εe)^{1 / 5}. Deck et al. (2013) adapted those two criteria to the case of two massive planets, finding little difference up to the numerical coefficients. However, they treat two different situations; the case of orbits initially circular and the case of two eccentric orbits. As in Mustill & Wyatt (2012), the eccentric criterion proposed can be used for eccentricities verifying e₁ + e₂ ≳ 1.33 ε^{3 / 7}. We show here that the two Deck’s criteria can be obtained as the limit cases of a general expression.

3.1. Width of the libration area

Using the same approach as Wisdom (1980), Deck et al. (2013), we have to express the width of the resonant island as a function of the orbital parameters and compare it with the distance between the two adjacent centers of MMR.

In the (X,Y) plane, the center of the resonance is located at the point of coordinates (X₁,0). The width of the libration area is defined as the distance between the two separatrices on the Y = 0 axis. It is indeed the direction where the resonant island is the widest.

We note ${\mathit{X}}_{\mathrm{1}}^{\mathrm{\ast }}\mathit{,}{\mathit{X}}_{\mathrm{2}}^{\mathrm{\ast }}$ the abscissas of the intersections between the separatrices and the Y = 0 axis. Relations between ${\mathit{X}}_{\mathrm{1}}^{\mathrm{\ast }}\mathit{,}{\mathit{X}}_{\mathrm{2}}^{\mathrm{\ast }}$, and X₃ can be derived (see Appendix C.1) and we obtain the expressions of ${\mathit{X}}_{\mathrm{1}}^{\mathrm{\ast }}$ and ${\mathit{X}}_{\mathrm{2}}^{\mathrm{\ast }}$ as functions of X₃ (Ferraz-Mello 2007; Deck et al. 2013). We have $\begin{array}{ccc}{\mathit{X}}_{\mathrm{1}}^{\mathrm{\ast }}& \mathrm{=}& \mathrm{-}{\mathit{X}}_{\mathrm{3}}\mathrm{-}\frac{\mathrm{2}}{\sqrt{{\mathit{X}}_{\mathrm{3}}}}\mathit{,}\\ {\mathit{X}}_{\mathrm{2}}^{\mathrm{\ast }}& \mathrm{=}& \mathrm{-}{\mathit{X}}_{\mathrm{3}}\mathrm{+}\frac{\mathrm{2}}{\sqrt{{\mathit{X}}_{\mathrm{3}}}}\mathrm{·}\end{array}$The width of the libration zone δX depends solely on the value of X₃, $\mathit{\delta X}\mathrm{=}\frac{\mathrm{4}}{\sqrt{{\mathit{X}}_{\mathrm{3}}}}\mathrm{·}$(55)In order to study the overlap of resonance islands, we need the width of the resonance in terms of α. Let us invert the previous change of variables in order to express the variation of α in terms of the variation of X. In this subsection, for any function Q(X), we note $\mathit{\delta Q}\mathrm{=}\mathrm{\left|}\mathit{Q}\mathrm{\right(}{\mathit{X}}_{\mathrm{1}}^{\mathrm{\ast }}\mathrm{)}\mathrm{-}\mathit{Q}\mathrm{(}{\mathit{X}}_{\mathrm{2}}^{\mathrm{\ast }}\mathrm{\left)}\mathrm{\right|}\mathit{.}$(56)The computation of δℐ (47) is straightforward from the computation of δX$\begin{array}{ccc}\mathit{\delta }\mathrm{ℐ}& \mathrm{=}& \left|\frac{{{\mathit{X}}_{\mathrm{1}}^{\mathrm{\ast }}}^{\mathrm{2}}}{\mathrm{2}}\mathrm{-}\frac{{{\mathit{X}}_{\mathrm{2}}^{\mathrm{\ast }}}^{\mathrm{2}}}{\mathrm{2}}\right|\\ & \mathrm{=}& \frac{\mathrm{1}}{\mathrm{2}}\left|{\mathit{X}}_{\mathrm{2}}^{\mathrm{\ast }}\mathrm{+}{\mathit{X}}_{\mathrm{1}}^{\mathrm{\ast }}\right|\left|{\mathit{X}}_{\mathrm{2}}^{\mathrm{\ast }}\mathrm{-}{\mathit{X}}_{\mathrm{1}}^{\mathrm{\ast }}\right|\\ & \mathrm{=}& {\mathit{X}}_{\mathrm{3}}\mathit{\delta X}\\ \mathit{\delta }\mathrm{ℐ}& \mathrm{=}& \mathrm{4}\sqrt{{\mathit{X}}_{\mathrm{3}}}\mathit{.}\end{array}$(57)We then directly deduce δI₁ = χ^{2 / 3}δℐ from Eq. (46). Since I₂ and ΔG are first integrals, the variation of Λ_i only depends on δI₁. And finally, since we have $\mathit{\alpha }\mathrm{=}{\left({\mathit{\gamma }}^{-1}\frac{{\mathrm{\Lambda }}_{\mathrm{1}}}{{\mathrm{\Lambda }}_{\mathrm{2}}}\right)}^{\mathrm{2}}\mathit{,}$(58)αcan be developed to the first order in (C−ΔG) thanks to Eq. (34). This development gives $\mathit{\alpha }\mathrm{=}{\mathit{\alpha }}_{\mathrm{0}}\left(\mathrm{1}\mathrm{-}\frac{\mathrm{2}\mathrm{(}{\mathit{\alpha }}_{\mathrm{0}}\mathrm{+}\mathit{\gamma }{\mathrm{)}}^{\mathrm{2}}}{\mathit{\gamma }{\mathit{\alpha }}_{\mathrm{0}}}\mathrm{(}\mathit{p}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{(}{\mathit{I}}_{\mathrm{1}}\mathrm{-}{\mathit{\chi }}^{\mathrm{2}\mathit{/}\mathrm{3}}{\mathrm{ℐ}}_{\mathrm{0}}\mathrm{)}\right)\mathit{.}$(59)The width of the resonance in terms of α is then directly related to X₃ through $\mathit{\delta \alpha }\mathrm{=}{\mathit{\alpha }}_{\mathrm{0}}\frac{\mathrm{8}{\mathit{r}}^{\mathrm{2}\mathit{/}\mathrm{3}}}{{\mathrm{3}}^{\mathrm{2}\mathit{/}\mathrm{3}}}{\mathit{\epsilon }}^{\mathrm{2}\mathit{/}\mathrm{3}}\mathrm{(}\mathit{p}\mathrm{+}\mathrm{1}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{3}}\sqrt{{\mathit{X}}_{\mathrm{3}}}\mathrm{+}\mathrm{o}\mathrm{\left(}{\mathit{\epsilon }}^{\mathrm{2}\mathit{/}\mathrm{3}}\mathrm{\right(}\mathit{p}\mathrm{+}\mathrm{1}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathrm{)}\mathit{.}$(60)The computation of the width of resonance is thus reduced to the computation of the root X₃ as a function of the parameters. It should also be remarked that at the first order, the width of resonance does not depend on the mass ratio γ.

3.2. Minimal AMD of a resonance

We are now interested in the overlap of adjacent resonant islands. Planets trapped in the chaotic zone created by the overlap will experience variations of their actions eventually leading to collisions.

For a configuration close to a given resonance p:p + 1, the AMD can evolve toward higher values if the original value places the system in a configuration above the inner separatrix, eventually leading the planets to collision or chaotic motion in case of MMR overlap. On the other hand, if the initial AMD of the planets forces them to remain in the inner circulation region of the overlapped MMR islands, the system will remain stable in regards to this criterion. Since C = I₁ + I₂, and I₂ is a first integral, we define the minimal AMD of a resonance¹C_min(p) as the minimal value of I₁ to enter the resonant island given ΔG−I₂. Two cases must be discussed:

• the point I₁ = 0 is already in the libration zone and then C_min = 0;

• the point I₁ = 0 is in the inner circulation zone and then we have

${\mathit{C}}_{\min }\mathrm{=}{\mathit{I}}_{\mathrm{1}}\mathrm{\left(}{\mathit{X}}_{\mathrm{2}}^{\mathrm{\ast }}\mathrm{\right)}\mathrm{=}\frac{{\mathit{\chi }}^{\mathrm{2}\mathit{/}\mathrm{3}}}{\mathrm{2}}{\left({\mathit{X}}_{\mathrm{3}}\mathrm{-}\frac{\mathrm{2}}{\sqrt{{\mathit{X}}_{\mathrm{3}}}}\right)}^{\mathrm{2}}\mathrm{·}$(61)In the second case, we have an implicit expression of X₃ depending on C_min${\mathit{\chi }}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{3}}\sqrt{\mathrm{2}{\mathit{C}}_{\min }}\mathrm{=}{\mathit{X}}_{\mathrm{3}}\mathrm{-}\frac{\mathrm{2}}{\sqrt{{\mathit{X}}_{\mathrm{3}}}}\mathit{,}$(62)where χ was defined in Eq. (40). In other words, there is a one-to-one correspondence between C_min Eq. (61) and the Hamiltonian parameter ℐ₀ for C_min> 0. The shape of the resonance island is completely described by C_min.

We can also use the definition of C_min to give an expression depending on the system parameters $\begin{array}{ccc}{\mathit{C}}_{\min }& \mathrm{=}& {\mathit{I}}_{\mathrm{1}}\mathrm{=}{\mathit{u}}_{\mathrm{1}}\mathit{u̅}\mathrm{1}\\ & \mathrm{=}& {\left|\frac{{\mathit{R}}_{\mathrm{1}}}{\mathit{R}}\sqrt{\frac{{\mathrm{\Lambda }}_{\mathrm{1}\mathit{,}\mathrm{0}}}{\mathrm{2}}}{X}_{\mathrm{1}}\mathrm{+}\frac{{\mathit{R}}_{\mathrm{2}}}{\mathit{R}}\sqrt{\frac{{\mathrm{\Lambda }}_{\mathrm{2}\mathit{,}\mathrm{0}}}{\mathrm{2}}}{X}_{\mathrm{2}}\right|}^{\mathrm{2}}\\ & \mathrm{=}& {\left(\frac{{\mathit{R}}_{\mathrm{1}}}{\mathit{R}}\right)}^{\mathrm{2}}\frac{{\mathrm{\Lambda }}_{\mathrm{1}\mathit{,}\mathrm{0}}}{\mathrm{2}}{\left|X_{\mathrm{1}}\mathrm{-}\left|\frac{{\mathit{R}}_{\mathrm{2}}}{{\mathit{R}}_{\mathrm{1}}}\right|\sqrt{\frac{{\mathrm{\Lambda }}_{\mathrm{2}\mathit{,}\mathrm{0}}}{{\mathrm{\Lambda }}_{\mathrm{1}\mathit{,}\mathrm{0}}}}{X}_{\mathrm{2}}\right|}^{\mathrm{2}}\\ & \mathrm{\simeq }& \end{array}$(63)where ${\mathit{c}}_{\mathit{i}}\mathrm{=}\sqrt{\mathrm{2}}\sqrt{\mathrm{1}\mathrm{-}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathit{i}}^{\mathrm{2}}}}\mathrm{=}\mathrm{\left|}{X}_{\mathit{i}}\mathrm{\right|}$. We note ${\mathit{c}}_{\min }\mathrm{=}{\mathit{c}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{\mathit{c}}_{\mathrm{2}}^{\mathrm{2}}\mathrm{-}\mathrm{2}{\mathit{c}}_{\mathrm{1}}{\mathit{c}}_{\mathrm{2}}\cos \mathrm{\Delta }\mathit{\varpi ,}$(64)the reduced minimal AMD. We can use the expression (63) to compute the quantity ${\mathit{\chi }}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{3}}\sqrt{\mathrm{2}{\mathit{C}}_{\min }}$ appearing in Eq. (62) ${\mathit{\chi }}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{3}}\sqrt{\mathrm{2}{\mathit{C}}_{\min }}\mathrm{\simeq }\frac{{\mathrm{3}}^{\mathrm{1}\mathit{/}\mathrm{3}}}{{\mathit{r}}^{\mathrm{1}\mathit{/}\mathrm{3}}}\frac{\mathrm{(}\mathit{p}\mathrm{+}\mathrm{1}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{3}}}{{\mathit{\epsilon }}^{\mathrm{1}\mathit{/}\mathrm{3}}}\sqrt{{\mathit{c}}_{\min }}\mathrm{+}\mathrm{o}\mathrm{\left(}{\mathit{p}}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathrm{\right)}\mathit{.}$(65)The function C_min(X₃) (Eq. (61)) is plotted in Fig. 2 with the two approximations used by Deck et al. (2013) to obtain the width of the resonance. For C_min ≫ χ^{2 / 3} or C_min close to zero, the relation can be simplified and we obtain $\begin{array}{ccc}{\mathit{X}}_{\mathrm{3}}& \mathrm{~}& \\ {\mathit{X}}_{\mathrm{3}}& \mathrm{=}& \end{array}$We can use the developments (66) and (67) in order to compute the width of the resonance in these two cases (see Appendix C). It should be noted as well that for C_min = 0, we have X₃ = 2^{2 / 3}.

Fig. 2 Relation Eq. (62) between X3 and Cmin Eq. (61) and two different approximations. In red, the approximation used by Deck et al. (2013) for eccentric orbits and in purple the constant evaluation used for circular orbits.

3.3. Implicit overlap criterion

The overlap of MMR can be determined by finding the first integer p such that the sum of the half-width of the resonances p:p + 1 and p + 1:p + 2 is larger than the distance between the respective centers of these two resonances (Wisdom 1980; Deck et al. 2013) $\frac{\mathrm{\Delta }\mathit{\alpha }}{{\mathit{\alpha }}_{\mathrm{0}\mathit{,p}}}\lesssim \frac{\mathrm{1}}{\mathrm{2}}\left(\frac{\mathit{\delta }{\mathit{\alpha }}_{\mathit{p}}}{{\mathit{\alpha }}_{\mathrm{0}\mathit{,p}}}\mathrm{+}\frac{\mathit{\delta }{\mathit{\alpha }}_{\mathit{p}\mathrm{+}\mathrm{1}}}{{\mathit{\alpha }}_{\mathrm{0}\mathit{,p}\mathrm{+}\mathrm{1}}}\right)\mathit{,}$(68)where Δα is the distance between the two centers and δα_k corresponds to the width of the resonance k:k + 1.

Up to terms of order ε^{2 / 3}, the center of the resonance island p:p + 1 is located at the center of the resonance of the unperturbed Keplerian problem, α_0,p = (p/ (p + 1))^{2 / 3}. We develop α_0,p for p ≫ 1$\begin{array}{ccc}{\mathit{\alpha }}_{\mathrm{0}\mathit{,p}}& \mathrm{=}& {\left(\frac{\mathit{p}}{\mathit{p}\mathrm{+}\mathrm{1}}\right)}^{\mathrm{2}\mathit{/}\mathrm{3}}\\ & \mathrm{=}& \end{array}$(69)Therefore, we have at second order in p$\frac{\mathrm{\Delta }\mathit{\alpha }}{{\mathit{\alpha }}_{\mathrm{0}\mathit{,p}}}\mathrm{=}\frac{\mathrm{2}}{\mathrm{3}}\frac{\mathrm{1}}{\mathrm{(}\mathit{p}\mathrm{+}\mathrm{1}{\mathrm{)}}^{\mathrm{2}}}\mathrm{·}$(70)We can use the implicit expression (62) of X₃ as a function of $\sqrt{{\mathit{c}}_{\min }}$ (Eq. (64)) in order to derive an overlap criterion independent of approximations on the value of C_min. Equating the general width of resonance (60) with the distance between to adjacent centers (70) and isolating X₃ gives ${\mathit{X}}_{\mathrm{3}}\mathrm{=}\frac{{\mathrm{3}}^{\mathrm{4}\mathit{/}\mathrm{3}}}{\mathrm{144}{\mathit{r}}^{\mathrm{4}\mathit{/}\mathrm{3}}}{\mathit{\epsilon }}^{\mathrm{-}\mathrm{4}\mathit{/}\mathrm{3}}\mathrm{(}\mathit{p}\mathrm{+}\mathrm{1}{\mathrm{)}}^{\mathrm{-}\mathrm{14}\mathit{/}\mathrm{3}}\mathit{.}$(71)We can inject this expression of X₃ into (62), and using Eq. (65), $\sqrt{{\mathit{c}}_{\min }}\mathrm{=}\frac{\mathrm{1}}{\mathrm{48}\mathit{r\epsilon }\mathrm{(}\mathit{p}\mathrm{+}\mathrm{1}{\mathrm{)}}^{\mathrm{5}}}\mathrm{-}\mathrm{8}\mathit{r\epsilon }\mathrm{(}\mathit{p}\mathrm{+}\mathrm{1}{\mathrm{)}}^{\mathrm{2}}\mathit{.}$(72)Using the first order expression of (p + 1) as a function of α, $\frac{\mathrm{1}}{\mathit{p}\mathrm{+}\mathrm{1}}\mathrm{=}\frac{\mathrm{3}}{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\alpha }\mathrm{)}$(73)we obtain an implicit expression of the overlap criterion $\sqrt{{\mathit{c}}_{\min }}\mathrm{=}\frac{{\mathrm{3}}^{\mathrm{4}}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\alpha }{\mathrm{)}}^{\mathrm{5}}}{{\mathrm{2}}^{\mathrm{9}}\mathit{r\epsilon }}\mathrm{-}\frac{\mathrm{32}\mathit{r\epsilon }}{\mathrm{9}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\alpha }{\mathrm{)}}^{\mathrm{2}}}\mathrm{·}$(74)

3.4. Overlap criterion for circular orbits

The implicit expression (74) can be used to find the criteria proposed by Deck et al. (2013) for circular and eccentric orbits. Let us first obtain the circular criterion by imposing c_min = 0 in Eq. (74) ${\mathrm{3}}^{\mathrm{6}}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\alpha }{\mathrm{)}}^{\mathrm{7}}\mathrm{=}{\mathrm{2}}^{\mathrm{14}}{\mathit{r}}^{\mathrm{2}}{\mathit{\epsilon }}^{\mathrm{2}}\mathit{.}$(75)We can express 1−α as a function of ε and we obtain $\mathrm{1}\mathrm{-}{\mathit{\alpha }}_{\mathrm{overlap}}\mathrm{=}\frac{\mathrm{4}{\mathit{r}}^{\mathrm{2}\mathit{/}\mathrm{7}}}{{\mathrm{3}}^{\mathrm{6}\mathit{/}\mathrm{7}}}{\mathit{\epsilon }}^{\mathrm{2}\mathit{/}\mathrm{7}}\mathrm{=}\mathrm{1.46}{\mathit{\epsilon }}^{\mathrm{2}\mathit{/}\mathrm{7}}\mathit{.}$(76)The exponent 2/7 was first proposed by Wisdom (1980) and the numerical factor 1.46 is similar to the one found by Deck et al. (2013).

Fig. 3 Representation of the MMR overlap criteria. The dotted lines correspond to the criteria proposed by Deck et al. (2013), and the collision curve is the approximation of the collision curve for α → 1. We represented in transparent green (p odd) and blue (p even) the first p:p + 1 MMR islands to show the agreement between the proposed overlap criterion and the actual intersections. In this figure, ε = 10-6.

3.5. Overlap criterion for high-eccentricity orbits

For large eccentricity, Deck et al. (2013) proposes a criterion based on the development (66) of Eq. (62). This criterion is obtained from Eq. (74) by ignoring the second term of the right-hand side which leads to ${\mathrm{2}}^{\mathrm{9}}\mathit{r\epsilon }\sqrt{{\mathit{c}}_{\min }}\mathrm{=}{\mathrm{3}}^{\mathrm{4}}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\alpha }{\mathrm{)}}^{\mathrm{5}}\mathit{.}$(77)Isolating 1−α gives $\mathrm{1}\mathrm{-}\mathit{\alpha }\mathrm{=}\frac{{\mathrm{2}}^{\mathrm{9}\mathit{/}\mathrm{5}}}{{\mathrm{3}}^{\mathrm{4}\mathit{/}\mathrm{5}}}{\mathit{r}}^{\mathrm{1}\mathit{/}\mathrm{5}}{\mathit{\epsilon }}^{\mathrm{1}\mathit{/}\mathrm{5}}{\mathit{c}}_{\min }^{\mathrm{1}\mathit{/}\mathrm{10}}\mathrm{=}\mathrm{1.38}{\mathit{\epsilon }}^{\mathrm{1}\mathit{/}\mathrm{5}}{\mathit{c}}_{\min }^{\mathrm{1}\mathit{/}\mathrm{10}}\mathit{.}$(78)This result is also similar to Deck’s one. For small c_min, the criterion (78) is less restrictive than the criterion (76) obtained for circular orbits. The comparison of these two overlap criteria provides a minimal value of c_min for the validity of the eccentric criterion $\sqrt{{\mathit{c}}_{\min }}\mathrm{=}\mathrm{1.33}{\mathit{\epsilon }}^{\mathrm{3}\mathit{/}\mathrm{7}}\mathit{.}$(79)

3.6. Overlap criterion for low-eccentricity orbits

For smaller eccentricities, we can develop Eq. (74) for small $\sqrt{{\mathit{c}}_{\min }}$ and α close to α_cir = 1−1.46ε^{2 / 7}, the critical semi major axis ratio for the circular overlap criterion (76). We have ${\mathrm{3}}^{\mathrm{2}}{\mathrm{2}}^{\mathrm{9}}\mathit{r\epsilon }\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\alpha }{\mathrm{)}}^{\mathrm{2}}\sqrt{{\mathit{c}}_{\min }}\mathrm{=}{\mathrm{3}}^{\mathrm{6}}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\alpha }{\mathrm{)}}^{\mathrm{7}}\mathrm{-}{\mathrm{2}}^{\mathrm{14}}{\mathit{r}}^{\mathrm{2}}{\mathit{\epsilon }}^{\mathrm{2}}\mathit{.}$(80)We develop the right-hand side at the first order in (α_cir−α) and evaluate the left-hand side for α = α_cir and after some simplifications obtain ${\mathit{\alpha }}_{\mathrm{cir}}\mathrm{-}\mathit{\alpha }\mathrm{=}\frac{{\mathrm{2}}^{\mathrm{9}}\mathit{r\epsilon }}{\mathrm{7}\mathrm{\times }{\mathrm{3}}^{\mathrm{4}}}\frac{\sqrt{{\mathit{c}}_{\min }}}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\alpha }}_{\mathrm{cir}}{\mathrm{)}}^{\mathrm{4}}}\mathrm{·}$(81)We inject the expression of α_cir into this equation and obtain the following development of the overlap criterion for low eccentricity: ${\mathit{\alpha }}_{\mathrm{cir}}\mathrm{-}\mathit{\alpha }\mathrm{=}\frac{\mathrm{2}\sqrt{{\mathit{c}}_{\min }}}{\mathrm{7}\mathrm{\times }{\mathrm{3}}^{\mathrm{4}\mathit{/}\mathrm{7}}{\mathit{r}}^{\mathrm{1}\mathit{/}\mathrm{7}}{\mathit{\epsilon }}^{\mathrm{1}\mathit{/}\mathrm{7}}}\mathrm{=}\mathrm{0.157}\frac{\sqrt{{\mathit{c}}_{\min }}}{{\mathit{\epsilon }}^{\mathrm{1}\mathit{/}\mathrm{7}}}\mathrm{·}$(82)This development remains valid for small enough $\sqrt{{\mathit{c}}_{\min }}$ if α_cir−α ≪ 1−α_cir, which can be rewritten $\mathrm{0.157}{\mathit{\epsilon }}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{7}}\sqrt{{\mathit{c}}_{\min }}\mathrm{\ll }\mathrm{1.46}{\mathit{\epsilon }}^{\mathrm{2}\mathit{/}\mathrm{7}}\mathit{,}$(83)which leads to $\sqrt{{\mathit{c}}_{\min }}\mathrm{\ll }\mathrm{9.30}{\mathit{\epsilon }}^{\mathrm{3}\mathit{/}\mathrm{7}}\mathit{.}$(84)It is worth noting that the low-eccentricity approximation allows to cover the range of eccentricities where the criterion (78) is not applicable, since both boundaries depend on the same power of ε.

We plot in Fig. 3 the overlap criteria (74) for ε = 10^-6, the two approximations (76) and (78) from Deck et al. (2013), as well as the collision condition used in Laskar & Petit (2017) approximated for α → 1, $\mathrm{1}\mathrm{-}\mathit{\alpha }\mathrm{\simeq }{\mathit{e}}_{\mathrm{1}}\mathrm{+}{\mathit{e}}_{\mathrm{2}}\mathrm{\simeq }\sqrt{{\mathit{c}}_{\min }}\mathit{.}$(85)We also plot the first MMR islands in order to show the agreement between the proposed criterion and the actual intersections. We see that for high eccentricities, and large 1−α, the system can verify the MMR overlap stability criterion while allowing for collision between the planets. For small α, the MMR overlap criterion alone cannot account for the stability of the system.

4. Critical AMD and MMR

4.1. Critical AMD in a context of resonance overlap

In Laskar & Petit (2017), we present the AMD-stability criterion based on the conservation of AMD. We assume the system dynamics to be secular chaotic. As a consequence the averaged semi-major axis and the total averaged AMD are conserved. Moreover, in this approximation the dynamics is limited to random AMD exchanges between planets with conservation of the total AMD. Based on these assumptions, collisions between planets are possible only if the AMD of the system can be distributed such that the eccentricities of the planets allow for collisions. Particularly, for each pair of adjacent planets, there exists a critical AMD, noted C_c(α,γ), such that for smaller AMD, collisions are forbidden.

The critical AMD was determined thanks to the limit collision condition $\mathit{\alpha }\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{e}}_{\mathrm{1}}\mathrm{)}\mathrm{=}\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{2}}\mathit{.}$(86)However, in practice, the system may become unstable long before orbit intersections; in particular the secular assumption does not hold if the system experiences chaos induced by MMR overlap. We can, though, consider that if the islands do not overlap, the AMD is, on average, conserved on timescales of order ε^{− 2 / 3} (i.e., of the order of the libration timescales). Therefore, the conservation, on average, of the AMD is ensured as long as the system adheres to the above criteria for any distribution of the AMD between planets. Based on the model of Laskar & Petit (2017), we compute a critical AMD associated to the criterion (74).

We consider a pair as AMD-stable if no distribution of AMD between the two planets allows the overlap of MMR. A first remark is that no pair can be considered as AMD-stable if α>α_cir, because in this case, even the circular orbits lead to MMR overlap. Let us write the criterion (74) as a function of α and ε; $\sqrt{{\mathit{c}}_{\min }}\mathrm{=}\mathit{g}\mathrm{\left(}\mathit{\alpha ,\epsilon }\mathrm{\right)}\mathit{,}$(87)where $\begin{array}{ccc}\mathit{g}\mathrm{\left(}\mathit{\alpha ,\epsilon }\mathrm{\right)}& \mathrm{=}& \frac{{\mathrm{3}}^{\mathrm{4}}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\alpha }{\mathrm{)}}^{\mathrm{5}}}{{\mathrm{2}}^{\mathrm{9}}\mathit{r\epsilon }}\mathrm{-}\frac{\mathrm{32}\mathit{r\epsilon }}{\mathrm{9}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\alpha }{\mathrm{)}}^{\mathrm{2}}}\mathit{\alpha }\mathit{<}{\mathit{\alpha }}_{\mathrm{cir}}\mathit{,}\\ & \mathrm{=}& \mathrm{0}\mathit{\alpha }\mathit{>}{\mathit{\alpha }}_{\mathrm{cir}}\mathit{.}\end{array}$(88)$\sqrt{{\mathit{c}}_{\min }}$depends on Δϖ and has a maximum for Δϖ = π. Since the variation of Δϖ does not affect the AMD of the system, we fix Δϖ = π since it is the least-favorable configuration. Therefore we have $\sqrt{{\mathit{c}}_{\min }}\mathrm{=}{\mathit{c}}_{\mathrm{1}}\mathrm{+}{\mathit{c}}_{\mathrm{2}}\mathit{.}$(89)We define the relative AMD of a pair of planets C and express it as a function of the variables c_i$C\mathrm{=}\frac{\mathit{C}}{{\mathrm{\Lambda }}_{\mathrm{2}}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}\mathit{\gamma }\sqrt{\mathit{\alpha }}{\mathit{c}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{{\mathit{c}}_{\mathrm{2}}^{\mathrm{2}}}^{\mathrm{)}}\mathit{.}$(90)The critical AMD ${\mathit{C}}_{\mathrm{c}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}$ associated to the overlap criterion (74) can be defined as the smallest value of relative AMD such that the conditions $\begin{array}{ccc}E\mathrm{\left(}{\mathit{c}}_{\mathrm{1}}\mathit{,}{\mathit{c}}_{\mathrm{2}}\mathrm{\right)}& \mathrm{=}& \\ C\mathrm{\left(}{\mathit{c}}_{\mathrm{1}}\mathit{,}{\mathit{c}}_{\mathrm{2}}\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}\mathit{\gamma }\sqrt{\mathit{\alpha }}{\mathit{c}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{{\mathit{c}}_{\mathrm{2}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{=}{\mathit{C}}_{\mathrm{c}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\end{array}$are verified by any couple (c₁,c₂). We represent this configuration in Fig. 4.

Fig. 4 MMR overlap criterion represented in the (c1,c2) plane.

As in Laskar & Petit (2017), the critical AMD is obtained through Lagrange multipliers $\mathrm{\nabla }C\mathrm{\propto }\mathrm{\nabla }E\mathit{.}$(93)The tangency condition gives a relation between c₁ and c₂, $\mathit{\gamma }\sqrt{\mathit{\alpha }}{\mathit{c}}_{\mathrm{1}}\mathrm{=}{\mathit{c}}_{\mathrm{2}}\mathit{.}$(94)Replacing c₂ in relation (91) gives the critical expression of c₁ and we immediately obtain the expression of c₂${\mathit{c}}_{\mathit{c,}\mathrm{1}}\mathrm{=}\frac{\mathit{g}\mathrm{\left(}\mathit{\alpha ,\epsilon }\mathrm{\right)}}{\mathrm{1}\mathrm{+}\mathit{\gamma }\sqrt{\mathit{\alpha }}} {\mathit{c}}_{\mathit{c,}\mathrm{2}}\mathrm{=}\frac{\mathit{\gamma }\sqrt{\mathit{\alpha }}\mathit{g}\mathrm{\left(}\mathit{\alpha ,\epsilon }\mathrm{\right)}}{\mathrm{1}\mathrm{+}\mathit{\gamma }\sqrt{\mathit{\alpha }}}\mathrm{·}$(95)The value of ${\mathit{C}}_{\mathrm{c}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}$ is obtained by injecting the critical values c_c,1 and c_c,2 into the expression of C${\mathit{C}}_{\mathrm{c}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\mathrm{\left(}\mathit{\alpha ,\gamma ,\epsilon }\mathrm{\right)}\mathrm{=}\frac{\mathit{g}\mathrm{(}\mathit{\alpha ,\epsilon }{\mathrm{)}}^{\mathrm{2}}}{\mathrm{2}}\frac{\mathit{\gamma }\sqrt{\mathit{\alpha }}}{\mathrm{1}\mathrm{+}\mathit{\gamma }\sqrt{\mathit{\alpha }}}\mathrm{·}$(96)

Fig. 5 Regions of application of the different criteria presented in this work. The purple region represents ${\mathit{C}}_{\mathrm{c}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}$ is the smallest, in the green zone, ${\mathit{C}}_{\mathrm{c}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}$ is the smallest and the circular overlap criterion is verified in the red zone. We see that the curve αR computed through a development of ${\mathit{C}}_{\mathrm{c}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}$ and ${\mathit{C}}_{\mathrm{c}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}$ presents a good agreement with the real limit between the green and the purple area. Here γ = 1.

4.2. Comparison with the collision criterion

It is then natural to compare the critical AMD ${\mathit{C}}_{\mathrm{c}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}$ to the critical AMD C_c (denoted hereafter by ${\mathit{C}}_{\mathrm{c}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}$) derived from the collision condition (86). If α>α_cir, the circular overlap criterion implies that ${\mathit{C}}_{\mathrm{c}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\mathrm{=}\mathrm{0}$ and therefore ${\mathit{C}}_{\mathrm{c}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}$ should be preferred to the previous criterion ${\mathit{C}}_{\mathrm{c}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}$. However, ${\mathit{C}}_{\mathrm{c}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}$ was obtained thanks to the assumption that α was close to 1. Particularly, it makes no sense to talk about first-order MMR overlap for α< 0.63 which corresponds to the center of the MMR 2:1. Therefore, the collision criterion should be used for small α. We need then to find α_R such that for α<α_R, we should use the critical AMD ${\mathit{C}}_{\mathrm{c}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}$. Since we are close to 1, we use a development of ${\mathit{C}}_{\mathrm{c}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}$ presented in Laskar & Petit (2017), and similarly, only keep the leading terms in 1−α in ${\mathit{C}}_{\mathrm{c}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}$. The two expressions are ${\mathit{C}}_{\mathrm{c}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{=}\frac{\mathit{\gamma }}{\mathrm{1}\mathrm{+}\mathit{\gamma }}\frac{\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\alpha }{\mathrm{)}}^{\mathrm{2}}}{\mathrm{2}}\mathit{,} {\mathit{C}}_{\mathrm{c}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\mathrm{=}\frac{\mathit{\gamma }}{\mathrm{1}\mathrm{+}\mathit{\gamma }}\frac{\mathit{g}\mathrm{(}\mathit{\alpha ,\epsilon }{\mathrm{)}}^{\mathrm{2}}}{\mathrm{2}}\mathrm{·}$(97)We observe that for α close to 1, the two expression have the same dependence on γ, therefore, α_R depends solely on ε. Simplifying ${\mathit{C}}_{\mathrm{c}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{=}{\mathit{C}}_{\mathrm{c}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}$ gives α_R as a solution of the polynomial equation in (1−α); ${\mathrm{3}}^{\mathrm{6}}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\alpha }{\mathrm{)}}^{\mathrm{7}}\mathrm{-}{\mathrm{3}}^{\mathrm{2}}{\mathrm{2}}^{\mathrm{9}}\mathit{r\epsilon }\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\alpha }{\mathrm{)}}^{\mathrm{3}}\mathrm{-}{\mathrm{2}}^{\mathrm{14}}\mathrm{(}\mathit{r\epsilon }{\mathrm{)}}^{\mathrm{2}}\mathrm{=}\mathrm{0.}$(98)While an exact analytical solution cannot be provided, a development in powers of ε gives the following expression $\begin{array}{ccc}\mathrm{1}\mathrm{-}{\mathit{\alpha }}_{\mathit{R}}& \mathrm{=}& \frac{\mathrm{4}}{\mathrm{3}}\mathrm{(}\mathrm{2}\mathit{r\epsilon }{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{4}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{4}}\sqrt{\mathrm{2}\mathit{r\epsilon }}\mathrm{+}\mathrm{O}\mathrm{(}{\mathit{\epsilon }}^{\mathrm{3}\mathit{/}\mathrm{4}}\mathrm{)}\\ & \mathrm{=}& \end{array}$(99)It should be remarked that the first term can be directly obtained using Deck’s high-eccentricity approximation.

In Fig. 5 we plot α_R and α_cir and indicate which criterion is used in the areas delimited by the curves. We specifically represented the region α>α_cir because we cannot treat this region in a similar manner to the remaining region since comparing the relative AMD C to ${\mathit{C}}_{\mathrm{c}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}$ does not provide any information. We see that the curve α_R is not exactly at the limit where ${\mathit{C}}_{\mathrm{c}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{=}{\mathit{C}}_{\mathrm{c}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}$ for higher ε due to the development of the critical AMDs for α → 1. We study the influence of γ on the difference between α_R and the actual limit in Appendix D

For stability analysis, we need to choose the smallest of the two critical AMD. For α<α_R, the collisional criterion is better and the MMR overlap criterion is used for α>α_R. We thus define a piece-wise global critical AMD represented in Fig. 6$\begin{array}{ccc}{\mathit{C}}_{\mathrm{c}}\mathrm{\left(}\mathit{\alpha ,\gamma ,\epsilon }\mathrm{\right)}& \mathrm{=}& {\mathit{C}}_{\mathrm{c}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{\left(}\mathit{\alpha ,\gamma }\mathrm{\right)}\mathit{\alpha }\mathit{<}{\mathit{\alpha }}_{\mathit{R}}\mathrm{\left(}\mathit{\epsilon ,\gamma }\mathrm{\right)}\mathit{,}\\ & \mathrm{=}& \end{array}$(100)

Fig. 6 Representation of the two critical AMD presented in this paper. ${\mathit{C}}_{\mathrm{c}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}$ in black is the collisional criterion from Laskar & Petit (2017), ${\mathit{C}}_{\mathrm{c}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}$ in red is the critical AMD derived from the MMR overlap criterion. In this plot, ε = 10-4 and γ = 1.

5. Effects of the MMR overlap on the AMD-classification of planetary systems

In Laskar & Petit (2017), we proposed a classification of the planetary systems based on their AMD-stability. A system is considered as AMD-stable if every adjacent pair of planets is AMD-stable. A pair is considered as AMD-stable if its AMD-stability coefficient $\mathit{\beta }\mathrm{=}\frac{\mathit{C}}{{\mathrm{\Lambda }}^{\mathrm{\prime }}{\mathit{C}}_{\mathrm{c}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}}\mathit{<}\mathrm{1}\mathit{,}$(101)where C is the total AMD of the system, Λ′ is the circular momentum of the outer planet and ${\mathit{C}}_{\mathrm{c}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}$ is the critical AMD derived from the collision condition. A similar AMD-coefficient can be defined using the global critical AMD defined in (100) instead of the collisional critical AMD ${\mathit{C}}_{\mathrm{c}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}$. Let us note β^(MMR), the AMD-stability coefficient associated to the critical AMD (100).

We can first observe that β^(MMR) is not defined for α>α_cir. Indeed, the conservation of the AMD cannot be guaranteed for orbits experiencing short-term chaos.

We use the modified definition of AMD-stability in order to test its effects on the AMD-classification proposed in Laskar & Petit (2017).

5.1. Sample and methodology

We first briefly recall the methodology used in Laskar & Petit (2017); to which we refer the reader for full details. We compute the AMD-stability coefficients for the systems taken from the Extrasolar Planets Encyclopaedia² with known periods, planet masses, eccentricities, and stellar mass. For each pair of adjacent planets, ε was computed using the expression $\mathit{\epsilon }\mathrm{=}\frac{{\mathit{m}}_{\mathrm{1}}\mathrm{+}{\mathit{m}}_{\mathrm{2}}}{{\mathit{m}}_{\mathrm{0}}}\mathit{,}$(102)where m₁ and m₂ are the two planet masses and m₀, the star mass. The semi-major axis ratio was derived from the period ratio and Kepler third law in order to reduce the uncertainty.

The systems are assumed coplanar, however in order to take into account the contribution of the real inclinations to the AMD, we define C_p, the coplanar AMD of the system, defined as the AMD of the same system if it was coplanar. We can compute coplanar AMD-stability coefficients ${\mathit{\beta }}_{\mathit{p}}^{\mathrm{\left(}\mathrm{MMR}\mathrm{\right)}}$ using C_p instead of C, and we define the total AMD-stability coefficients as $\mathit{\beta }\mathrm{=}\mathrm{2}{\mathit{\beta }}_{\mathit{p}}^{\mathrm{\left(}\mathrm{MMR}\mathrm{\right)}}$. Doing so, we assume the equipartition of the AMD between the different degree of freedom of the system.

We assume the uncertainties of the database quantities to be Gaussian. For the eccentricities, we use the same method as in the previous paper. The quantity ecosϖ is assumed to be Gaussian with the mean, the value of the database and standard deviation, the database uncertainty. The quantity esinϖ is assumed to have a Gaussian distribution with zero mean and the same standard deviation. The distribution of eccentricity is then derived from these two distributions.

We then propagate the uncertainties through the computations thanks to Monte Carlo simulations of the original distributions. For each of the systems, we drew 10 000 values of masses, periods and eccentricities from the computed distributions. We then compute β^(MMR) for each of these configurations and compute the 1-σ confidence interval.

In Laskar & Petit (2017), we studied 131 systems but we did not find the stellar mass for 4 of these systems. They were, as a consequence, excluded from this study. Moreover, the computation of ε for the pairs of planets of the 127 remaining systems of the sample led in some cases to high planet-to-star mass ratios. We decide to exclude the systems such that α_cir was smaller than the center of the resonance 2:1. We thus discard systems such that a pair of planets has $\mathit{\epsilon }\mathit{>}{\mathit{\epsilon }}_{\lim }\mathrm{=}\mathrm{8.20}\mathrm{\times }{\mathrm{10}}^{-3}\mathit{.}$(103)As a result, we only consider in this study 111 systems that meet the above requirements.

A pair is considered stable if the 1-σ confidence interval (84% of the simulated systems) of the AMD-stability coefficient β^(MMR) is below 1. A system is stable if all adjacent pairs are stable.

5.2. Results