© ESO, 2017

1. Introduction

The increasing number of planetary systems has made it necessary to search for a possible classification of these planetary systems. Ideally, such a classification should not require heavy numerical analysis as it needs to be applied to large sets of systems. Some possible approaches can rely on the stability analysis of these systems, as this stability analysis is also part of the process used to consolidate the discovery of planetary systems. The stability analysis can also be considered a key part to understanding the wider question of the architecture of planetary systems. In fact, the distances between planets and other orbital characteristic distributions is one of the oldest questions in celestial mechanics, the most famous attempts to set laws for this distribution of planetary orbits being the Titius-Bode power laws (see Nieto 1972; Graner & Dubrulle 1994 for a review).

Recent research has focused on statistical analysis of observed architecture (Fabrycky et al. 2014; Lissauer et al. 2011; Mayor et al. 2011), eccentricity distribution (Moorhead et al. 2011; Shabram et al. 2016; Xie et al. 2016), or inclination distribution (Fang & Margot 2012; Figueira et al. 2012); see Winn & Fabrycky (2015) for a review. The analysis of these observations has been compared with models of system architecture (Fang & Margot 2013; Pu & Wu 2015; Tremaine 2015). These theoretical works usually use empirical criteria based on the Hill radius proposed by Gladman (1993) and refined by Chambers et al. (1996), Smith & Lissauer (2009), and Pu & Wu (2015). These criteria of stability usually multiply the Hill radius by a numerical factor Δ_sep empirically evaluated to a value around 10. They are extensions of the analytical results on Hill spheres for the three-body problem (Marchal & Bozis 1982). Works on chaotic motion caused by the overlap of mean motion resonances (MMR, Wisdom 1980; Deck et al. 2013; Ramos et al. 2015) could justify the Hill-type criteria, but the results on the overlap of the MMR island are valid only for close orbits and for short-term stability.

Another approach to stability analysis is to consider the secular approximation of a planetary system. In this framework, the conservation of the semi-major axis leads to the conservation of another quantity, the angular momentum deficit (AMD; Laskar 1997, 2000). An architecture model can be developed from this consideration (Laskar 2000). The AMD can be interpreted as a measure of the orbits’ excitation (Laskar 1997) that limits the planets close encounters and ensures long-term stability. Therefore a stability criterion can be derived from the semi-major axis, the masses and the AMD of a system. In addition, it can be demonstrated that the AMD decreases during inelastic collisions (see Sect. 2.3), accounting for the gain of stability of a lower multiplicity system. Here we extend the previous analysis of Laskar (2000), and derive more precisely the AMD-stability criterion that can be used to establish a classification of the multiplanetary systems.

In the original letter (Laskar 2000), the detailed computations were referred to as a preprint to be published. Although this preprint has been in nearly final form for more than a decade, and has even been provided to some researchers (Hernández-Mena & Benet 2011), it is still unpublished. In Sects. 2 and 3 we provide the fundamental concepts of AMD, the full description, and all proofs for the model that was described in Laskar (2000). This material is close to the unpublished preprint. Section 4 recalls briefly the model of planetary accretion based on AMD stochastic transfers that was first presented in Laskar (2000). This model provides analytical expressions for the averaged systems architecture and orbital parameter distribution, depending on the initial mass distribution (Table 2).

In Sect. 5, we show how the AMD-stability criterion presented in Sect. 3 can be used to develop a classification of planetary systems. This AMD-stability classification is then applied to a selection of 131 multiplanet systems from The Extrasolar Planet Encyclopaedia database (exoplanet.eu) with known eccentricities. Finally, in Sect. 6 we discuss our results and provide some possible extension for this work.

2. Angular momentum deficit

2.1. Planetary Hamiltonian

Let P₀,P₁,...,P_n be n + 1 bodies of masses m₀,m₁,...,m_n in gravitational interaction, and let O be their centre of mass. For every body P_i, we denote u_i = OP_i. In the barycentric reference frame with origin O, Newton’s equations of motion form a differential system of order 6(n + 1) and can be written in Hamiltonian form using the canonical coordinates $\mathrm{(}{{u}}_{\mathit{i}}\mathit{,}{\begin{array}{c}{}_{{˜}}\\ {u}\end{array}}_{\mathit{i}}\mathrm{=}{\mathit{m}}_{\mathit{i}}\dot{{u}}{}_{\mathit{i}}\mathrm{)}_{\mathit{i}\mathrm{=}\mathrm{0}\mathit{,n}}$ with Hamiltonian $\mathit{H}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\sum _{\mathit{i}\mathrm{=}\mathrm{0}}^{\mathit{n}}\frac{\mathrm{\parallel }{\begin{array}{c}{˜}\\ {u}\end{array}}_{\mathit{i}}{\mathrm{\parallel }}^{\mathrm{2}}}{{\mathit{m}}_{\mathit{i}}}\mathrm{-}\mathrm{𝒢}\sum _{\mathrm{0}\mathrm{\le }\mathit{i}\mathit{<}\mathit{j}}\frac{{\mathit{m}}_{\mathit{i}}{\mathit{m}}_{\mathit{j}}}{{\mathrm{\Delta }}_{\mathit{ij}}}\mathit{,}$(1)where Δ_ij = ∥ u_i − u_j ∥, and is the constant of gravitation. The reduction of the centre of mass is achieved by using the canonical heliocentric variables of Poincaré $\mathrm{\left(}{{r}}_{\mathit{i}}\mathit{,}{\begin{array}{c}{}_{{˜}}\\ {r}\end{array}}_{\mathit{i}}\mathrm{\right)}$ (Poincaré 1905; Laskar & Robutel 1995), defined as $\begin{array}{ccc}& & {{r}}_{\mathrm{0}}\mathrm{=}{{u}}_{\mathrm{0}}\mathrm{;} {{r}}_{\mathit{i}}\mathrm{=}{{u}}_{\mathit{i}}\mathrm{-}{{u}}_{\mathrm{0}} \mathrm{for}\mathit{i}\ne \mathrm{0}\\ & & {\begin{array}{c}{˜}\\ {r}\end{array}}_{\mathrm{0}}\mathrm{=}{\begin{array}{c}{˜}\\ {u}\end{array}}_{\mathrm{0}}\mathrm{+}{\begin{array}{c}{˜}\\ {u}\end{array}}_{\mathrm{1}}\mathrm{+}\mathrm{·}\mathrm{·}\mathrm{·}\mathrm{+}{\begin{array}{c}{˜}\\ {u}\end{array}}_{\mathit{n}}\mathrm{;} {\begin{array}{c}{˜}\\ {r}\end{array}}_{\mathit{i}}\mathrm{=}{\begin{array}{c}{˜}\\ {u}\end{array}}_{\mathit{i}} \mathrm{for}\mathit{i}\ne \mathrm{0.}\end{array}$This Hamiltonian can then be split into an integrable part, H₀, and a perturbation, H₁, $\mathit{H}\mathrm{=}{\mathit{H}}_{\mathrm{0}}\mathrm{+}\mathit{H}1$(4)with ${\mathit{H}}_{\mathrm{0}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}\frac{\mathrm{\parallel }{\begin{array}{c}{˜}\\ {r}\end{array}}_{\mathit{i}}{\mathrm{\parallel }}^{\mathrm{2}}}{{\mathit{m}}_{\mathit{i}}}\mathrm{-}\mathrm{𝒢}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}\frac{{\mathit{m}}_{\mathrm{0}}{\mathit{m}}_{\mathit{i}}}{{\mathit{r}}_{\mathit{i}}}$(5)and ${\mathit{H}}_{\mathrm{1}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\frac{\mathrm{\parallel }{\begin{array}{c}{˜}\\ {u}\end{array}}_{\mathrm{0}}{\mathrm{\parallel }}^{\mathrm{2}}}{{\mathit{m}}_{\mathrm{0}}}\mathrm{-}\mathrm{𝒢}\sum _{\mathrm{1}\mathrm{\le }\mathit{i}\mathit{<}\mathit{j}}\frac{{\mathit{m}}_{\mathit{i}}{\mathit{m}}_{\mathit{j}}}{{\mathrm{\Delta }}_{\mathit{ij}}}\mathrm{·}$(6)The integrable part, H₀, is the Hamiltonian of a sum of disjoint Kepler problems of a single planet of mass m_i around a fixed Sun of mass m₀. A set of adapted variables for H₀ will thus be given by the elliptical elements, (a_k,e_k,i_k,λ_k,ϖ_k,Ω_k), where a_k is the semi-major axis, e_k the eccentricity, i_k the inclination, λ_k the mean longitude, ϖ_k the longitude of the perihelion, and Ω_k the longitude of the node. They are defined as the elliptical elements associated to the Hamiltonian ${\mathit{H}}_{\mathrm{0}\mathit{k}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\frac{\mathrm{\parallel }{\begin{array}{c}{˜}\\ {r}\end{array}}_{\mathit{k}}{\mathrm{\parallel }}^{\mathrm{2}}}{{\mathit{m}}_{\mathit{k}}}\mathrm{-}\mathit{\mu }\frac{{\mathit{m}}_{\mathit{k}}}{{\mathit{r}}_{\mathit{k}}}$(7)with .

2.2. Angular momentum deficit (AMD)

The total linear momentum of the system is null in the barycentric reference frame $\mathit{L}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{0}}^{\mathit{n}}{\mathit{m}}_{\mathit{i}}\dot{{u}}\mathit{i}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{0}}^{\mathit{n}}{\begin{array}{c}{˜}\\ {u}\end{array}}_{\mathit{i}}\mathrm{=}{\begin{array}{c}{˜}\\ {r}\end{array}}_{\mathrm{0}}\mathrm{=}\mathrm{0.}$(8)Let G be the total angular momentum. Its expression is not changed by the linear symplectic change of variable $\mathrm{\left(}{u}\mathit{,}\begin{array}{c}{}_{{˜}}\\ {u}\end{array}\mathrm{\right)}\mathrm{-}\mathrm{\to }\mathrm{\left(}{r}\mathit{,}\begin{array}{c}{}_{{˜}}\\ {r}\end{array}\mathrm{\right)}$. We have thus ${G}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{0}}^{\mathit{n}}{{u}}_{\mathit{i}}\hspace{0.17em}\mathrm{\wedge }{\begin{array}{c}{˜}\\ {u}\end{array}}_{\mathit{i}}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{{r}}_{\mathit{i}}\hspace{0.17em}\mathrm{\wedge }{\begin{array}{c}{˜}\\ {r}\end{array}}_{\mathit{i}}\mathit{.}$(9)When expressed in heliocentric variables, the angular momentum is thus the sum of the angular momentum of the Keplerian problems of the unperturbed Hamiltonian H₀. In particular, if the angular momentum direction is assumed to be the axis z, the norm of the angular momentum is $\mathit{G}\mathrm{=}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{n}}{\mathrm{\Lambda }}_{\mathit{k}}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathit{k}}^{\mathrm{2}}}\cos {\mathit{i}}_{\mathit{k}}\mathit{,}$(10)where ${\mathrm{\Lambda }}_{\mathit{k}}\mathrm{=}{\mathit{m}}_{\mathit{k}}\sqrt{\mathit{\mu }\hspace{0.17em}{\mathit{a}}_{\mathit{k}}}$. For such a system, the AMD is defined as the difference between the norm of the angular momentum of a coplanar and circular system with the same semi-major axis values and the norm of the angular momentum (G), i.e. (Laskar 1997, 2000) $\mathit{C}\mathrm{=}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{n}}{\mathrm{\Lambda }}_{\mathit{k}}\left(\mathrm{1}\mathrm{-}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathit{k}}^{\mathrm{2}}}\cos {\mathit{i}}_{\mathit{k}}\right)\mathit{.}$(11)

2.3. AMD and collision of orbits

The instabilities of a planetary system often result in a modification of its architecture. A planet can be ejected from the system or can fall into the star; in both cases this results in a loss of AMD for the system. The outcome of the AMD after a planetary collisions is less trivial and needs to be computed. Assume that among our n + 1 bodies, the (totally inelastic) collision of two bodies of masses m₁ and m₂, and orbits occurs, forming a new body . During this collision we assume that the other bodies are not affected. The mass is conserved ${\mathit{m}}_{\mathrm{3}}\mathrm{=}{\mathit{m}}_{\mathrm{1}}\mathrm{+}{\mathit{m}}_{\mathrm{2}}\mathit{,}$(12)and the linear momentum in the barycentric reference frame is conserved so ${\begin{array}{c}{}_{{˜}}\\ {u}\end{array}}_{\mathrm{3}}\mathrm{=}{\begin{array}{c}{}_{{˜}}\\ {u}\end{array}}_{\mathrm{1}}\mathrm{+}{\begin{array}{c}{}_{{˜}}\\ {u}\end{array}}_{\mathrm{2}}$, and also ${\begin{array}{c}{˜}\\ {r}\end{array}}_{\mathrm{3}}\mathrm{=}{\begin{array}{c}{˜}\\ {r}\end{array}}_{\mathrm{1}}\mathrm{+}{\begin{array}{c}{˜}\\ {r}\end{array}}_{\mathrm{2}}\mathit{.}$(13)On the other hand, r₃ = r₁ = r₂ at the time of the collision, so the angular momentum is also conserved ${{r}}_{\mathrm{3}}\hspace{0.17em}\mathrm{\wedge }{\begin{array}{c}{˜}\\ {r}\end{array}}_{\mathrm{3}}\mathrm{=}{{r}}_{\mathrm{1}}\hspace{0.17em}\mathrm{\wedge }{\begin{array}{c}{˜}\\ {r}\end{array}}_{\mathrm{1}}\mathrm{+}{{r}}_{\mathrm{2}}\hspace{0.17em}\mathrm{\wedge }{\begin{array}{c}{˜}\\ {r}\end{array}}_{\mathrm{2}}\mathrm{·}$(14)It should be noted that the transformation of the orbits during the collision is perfectly defined by Eqs. (12), (13). The problem which remains is to compute the evolution of the elliptical elements during the collision.

2.3.1. Energy evolution during collision

Just before the collision, we assume that the orbits and are elliptical heliocentric orbits. At the time of the collision, only these two bodies are involved and the other bodies are not affected. The evolution of the orbits are thus given by the conservation laws Eqs. (12), (13). The Keplerian energy of each particle is ${\mathit{h}}_{\mathit{i}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\frac{\mathrm{\parallel }{\begin{array}{c}{˜}\\ {r}\end{array}}_{\mathit{i}}{\mathrm{\parallel }}^{\mathrm{2}}}{{\mathit{m}}_{\mathit{i}}}\mathrm{-}\mathit{\mu }\frac{{\mathit{m}}_{\mathit{i}}}{{\mathit{r}}_{\mathit{i}}}\mathrm{=}\mathrm{-}\mathit{\mu }\frac{{\mathit{m}}_{\mathit{i}}}{\mathrm{2}{\mathit{a}}_{\mathit{i}}}\mathrm{·}$(15)At collision, we have r₁ = r₂ = r₃ = r; we have thus the conservation of the potential energy $\mathrm{-}\mathit{\mu }\frac{{\mathit{m}}_{\mathrm{3}}}{{\mathit{r}}_{\mathrm{3}}}\mathrm{=}\mathrm{-}\mathit{\mu }\frac{{\mathit{m}}_{\mathrm{1}}\mathrm{+}{\mathit{m}}_{\mathrm{2}}}{\mathit{r}}\mathrm{=}\mathrm{-}\mathit{\mu }\frac{{\mathit{m}}_{\mathrm{1}}}{{\mathit{r}}_{\mathrm{1}}}\mathrm{-}\mathit{\mu }\frac{{\mathit{m}}_{\mathrm{2}}}{{\mathit{r}}_{\mathrm{2}}}\mathrm{·}$(16)The change of Keplerian energy is thus given by the change of kinetic energy $\mathit{\delta h}\mathrm{=}{\mathit{h}}_{\mathrm{3}}\mathrm{-}{\mathit{h}}_{\mathrm{1}}\mathrm{-}{\mathit{h}}_{\mathrm{2}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\frac{\mathrm{\parallel }{\begin{array}{c}{˜}\\ {r}\end{array}}_{\mathrm{3}}{\mathrm{\parallel }}^{\mathrm{2}}}{{\mathit{m}}_{\mathrm{3}}}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\frac{\mathrm{\parallel }{\begin{array}{c}{˜}\\ {r}\end{array}}_{\mathrm{1}}{\mathrm{\parallel }}^{\mathrm{2}}}{{\mathit{m}}_{\mathrm{1}}}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\frac{\mathrm{\parallel }{\begin{array}{c}{˜}\\ {r}\end{array}}_{\mathrm{2}}{\mathrm{\parallel }}^{\mathrm{2}}}{{\mathit{m}}_{\mathrm{2}}}\mathrm{;}$(17)that is, with Eqs. (12), (13), $\begin{array}{ccc}& \mathrm{2}{\mathit{m}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}{\mathit{m}}_{\mathrm{3}}\mathit{\delta h}& \\ & \mathrm{=}{\mathit{m}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}\mathrm{(}{\begin{array}{c}{˜}\\ {r}\end{array}}_{\mathrm{1}}\mathrm{+}{\begin{array}{c}{˜}\\ {r}\end{array}}_{\mathrm{2}}{\mathrm{)}}^{\mathrm{2}}\mathrm{-}{\mathit{m}}_{\mathrm{2}}\mathrm{(}{\mathit{m}}_{\mathrm{1}}\mathrm{+}{\mathit{m}}_{\mathrm{2}}\mathrm{)}{\begin{array}{c}{˜}\\ {r}\end{array}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{-}{\mathit{m}}_{\mathrm{1}}\mathrm{(}{\mathit{m}}_{\mathrm{1}}\mathrm{+}{\mathit{m}}_{\mathrm{2}}\mathrm{)}{\begin{array}{c}{˜}\\ {r}\end{array}}_{\mathrm{2}}^{\mathrm{2}}& \\ & \mathrm{=}\mathrm{-}{\mathit{m}}_{\mathrm{2}}^{\mathrm{2}}{\begin{array}{c}{˜}\\ {r}\end{array}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{-}{\mathit{m}}_{\mathrm{1}}^{\mathrm{2}}{\begin{array}{c}{˜}\\ {r}\end{array}}_{\mathrm{2}}^{\mathrm{2}}\mathrm{+}\mathrm{2}{\mathit{m}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}{\begin{array}{c}{˜}\\ {r}\end{array}}_{\mathrm{1}}\mathrm{·}{\begin{array}{c}{˜}\\ {r}\end{array}}_{\mathrm{2}}& \\ & \mathrm{=}\mathrm{-}\mathrm{(}{\mathit{m}}_{\mathrm{2}}{\begin{array}{c}{˜}\\ {r}\end{array}}_{\mathrm{1}}\mathrm{-}{\mathit{m}}_{\mathrm{1}}{\begin{array}{c}{˜}\\ {r}\end{array}}_{\mathrm{2}}{\mathrm{)}}^{\mathrm{2}}& \\ & \mathrm{\le }\mathrm{0.}& \end{array}$(18)Thus, the Keplerian energy of the system decreases during collision. Part of the kinetic energy is dispelled during the collision. As expected, there is no loss of energy when ${\mathit{m}}_{\mathrm{2}}{\begin{array}{c}{}_{{˜}}\\ {r}\end{array}}_{\mathrm{1}}\mathrm{=}{\mathit{m}}_{\mathrm{1}}{\begin{array}{c}{}_{{˜}}\\ {r}\end{array}}_{\mathrm{2}}$, that is, as m₁m₂ ≠ 0, when $\dot{{u}}{}_{\mathrm{1}}\mathrm{=}\dot{{u}}{}_{\mathrm{2}}$. As an immediate consequence of the decrease of energy during the collision, we have $\frac{\mathrm{1}}{{\mathit{a}}_{\mathrm{3}}}\mathrm{\ge }\frac{\mathit{m}}{{\mathit{a}}_{\mathrm{1}}}\mathrm{+}\frac{\mathrm{1}\mathrm{-}\mathit{m}}{{\mathit{a}}_{\mathrm{2}}}$(19)with $\mathit{m}\mathrm{=}\frac{{\mathit{m}}_{\mathrm{1}}}{{\mathit{m}}_{\mathrm{1}}\mathrm{+}{\mathit{m}}_{\mathrm{2}}}\mathrm{=}\frac{{\mathit{m}}_{\mathrm{1}}}{{\mathit{m}}_{\mathrm{3}}}\mathrm{·}$(20)

2.3.2. AMD evolution during collision

Let $\mathit{f}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\mathrm{1}\mathit{/}\sqrt{\mathit{x}}$. As f′(x) < 0 and f′′(x) > 0, we have, as f is decreasing and convex, $\mathit{f}\left(\frac{\mathrm{1}}{{\mathit{a}}_{\mathrm{3}}}\right)\mathrm{\le }\mathit{f}\left(\frac{\mathit{m}}{{\mathit{a}}_{\mathrm{1}}}\mathrm{+}\frac{\mathrm{1}\mathrm{-}\mathit{m}}{{\mathit{a}}_{\mathrm{2}}}\right)\mathrm{\le }\mathit{mf}\left(\frac{\mathrm{1}}{{\mathit{a}}_{\mathrm{1}}}\right)\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{m}\mathrm{)}\mathit{f}\left(\frac{\mathrm{1}}{{\mathit{a}}_{\mathrm{2}}}\right)$(21)thus ${\mathit{m}}_{\mathrm{3}}\sqrt{{\mathit{a}}_{\mathrm{3}}}\mathrm{\le }{\mathit{m}}_{\mathrm{1}}\sqrt{{\mathit{a}}_{\mathrm{1}}}\mathrm{+}{\mathit{m}}_{\mathrm{2}}\sqrt{{\mathit{a}}_{\mathrm{2}}}\mathit{.}$(22)During the collision, the angular momentum is conserved Eq. (14), and so is the conservation of its normal component, that is $\begin{array}{ccc}& & {\mathit{m}}_{\mathrm{3}}\sqrt{\mathit{\mu }{\mathit{a}}_{\mathrm{3}}}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{3}}^{\mathrm{2}}}\cos {\mathit{i}}_{\mathrm{3}}\mathrm{=}\\ & & \end{array}$(23)We deduce that in all circumstances we have a decrease of the angular momentum deficit during the collision, that is ${\mathit{C}}_{\mathrm{3}}\mathrm{\le }{\mathit{C}}_{\mathrm{1}}\mathrm{+}{\mathit{C}}_{\mathrm{2}}$(24)with ${\mathit{C}}_{\mathit{k}}\mathrm{=}{\mathit{m}}_{\mathit{k}}\sqrt{\mathit{\mu }{\mathit{a}}_{\mathit{k}}}\mathrm{(}\mathrm{1}\mathrm{-}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathit{k}}^{\mathrm{2}}}\cos {\mathit{i}}_{\mathit{k}}\mathrm{)} \mathrm{(}\mathit{k}\mathrm{=}\mathrm{1}\mathit{,}\mathrm{3}\mathrm{)}\mathit{.}$(25)The equality can hold in Eq. (24) only if m₁ = 0, m₂ = 0, or a₁ = a₂ and , that is when one of the bodies is massless, or when the two bodies are on the same orbit and at the same position (at the time of the collision, we also have r₁ = r₂).

The diminution of AMD during collisions acts as a stabilisation of the system. A parallel can be made with thermodynamics, the AMD behaving for the orbits like the kinetic energy for the molecules of a perfect gas. The loss of AMD during collisions can thus be interpreted as a cooling of the system.

3. AMD-stability

We say that a planetary system is AMD-stable if the angular momentum deficit (AMD) amount in the system is not sufficient to allow for planetary collisions. As this quantity is conserved in the secular system at all orders (see Appendix B), we conjecture that in absence of short period resonances, the AMD-stability ensures the practical¹ long-term stability of the system. Thus for an AMD-stable system, short-term stability will imply long-term stability.

The condition of AMD-stability is obtained when the orbits of two planets of semi-major axis a,a′ cannot intersect under the assumption that the total AMD C has been absorbed by these two planets. It can be seen easily that the limit condition of collision is obtained in the planar case and can thus be written as $\begin{array}{ccc}& & \\ & & \end{array}$where (m,a,e) are the parameters of the inner orbit and (m′,a′,e′) those of the outer orbit (a ≤ a′).

3.1. Collisional condition

We assume that a,a′,m,m′ are non-zero. Denoting α = a/a′, ≳ = m/m′, the system in Eq. (27) becomes $\begin{array}{ccc}\mathrm{𝒟}\mathrm{\left(}\mathit{e,}{\mathit{e}}^{\mathrm{\prime }}\mathrm{\right)}& \mathrm{=}& \end{array}$(28)$\begin{array}{ccc}C\mathrm{\left(}\mathit{e,}{\mathit{e}}^{\mathrm{\prime }}\mathrm{\right)}& \mathrm{=}& \mathit{\gamma }\sqrt{\mathit{\alpha }}\mathrm{(}\mathrm{1}\mathrm{-}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}}\mathrm{)}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{\prime }\mathrm{2}}}\mathrm{)}\\ & \mathrm{=}& \end{array}$(29)with ${\mathrm{\Lambda }}^{\mathrm{\prime }}\mathrm{=}{\mathit{m}}^{\mathrm{\prime }}\sqrt{\mathit{\mu }{\mathit{a}}^{\mathrm{\prime }}}$, and where C(e,e′) = C/ Λ′ is called the relative AMD. As e and e′ are planetary eccentricities, we also have $\mathrm{0}\mathrm{\le }\mathit{e}\mathrm{\le }\mathrm{1}\mathrm{;} \mathrm{0}\mathrm{\le }{\mathit{e}}^{\mathrm{\prime }}\mathrm{\le }\mathrm{1.}$(30)

Fig. 1 Collision conditions for e0 = 1 /α − 1. Case (a): α< 1 / 2 ⇐ ⇒ e0> 1. Case (b): α> 1 / 2 ⇐ ⇒ e0< 1. Collisions can only occur in the shaded region.

The set of collisional conditions (Eqs. (28)–(30)) can be solved using Lagrange multipliers. We are looking for the minimum value of the relative AMD C(e,e′) (29) for which the collision condition (28) is satisfied. These conditions are visualised in the (e,e′) plane in Fig. 1. The collision condition (28) corresponds to the segment AB of Fig. 1. The domain of collisions is limited by the conditions (30). For e = 0, we have ${\mathit{e}}_{\mathrm{0}}^{\mathrm{\prime }}\mathrm{=}\mathrm{1}\mathrm{-}\mathit{\alpha }\mathit{<}\mathrm{1}$, and the intersection of the collision line with the axis e′ = 0 is obtained for e₀ = 1 /α − 1. This value can be greater or smaller than 1 depending on the value of α. We have thus the different cases (Fig. 1) $\begin{array}{ccc}\mathrm{\left(}\mathit{a}\mathrm{\right)} \mathit{\alpha }& \mathit{<}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\Leftarrow }\mathrm{\Rightarrow }{\mathit{e}}_{\mathrm{0}}\mathit{>}\mathrm{1} & \mathrm{and} \mathit{e}\mathrm{\le }\mathrm{1}\\ \mathrm{\left(}\mathit{b}\mathrm{\right)} \mathit{\alpha }& \mathit{>}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\Leftarrow }\mathrm{\Rightarrow }{\mathit{e}}_{\mathrm{0}}\mathit{<}\mathrm{1} & \mathrm{and} \mathit{e}\mathit{<}{\mathit{e}}_{\mathrm{0}}\end{array}$and the limit case α = 1 / 2, for which e₀ = 1. In all cases, the Lagrange multipliers condition is written $\mathit{\lambda }\mathrm{\nabla 𝒟}\mathrm{\left(}\mathit{e,}{\mathit{e}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{=}\mathrm{\nabla }C\mathrm{\left(}\mathit{e,}{\mathit{e}}^{\mathrm{\prime }}\mathrm{\right)}\mathit{,}$(33)which gives $\frac{\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{\prime }\mathrm{2}}}}{{\mathit{e}}^{\mathrm{\prime }}}\mathrm{=}\frac{\sqrt{\mathit{\alpha }}}{\gtrsim }\frac{\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}}}{\mathit{e}}\mathrm{·}$(34)This relation allows e′ to be eliminated in the collision condition (28), which becomes an equation in the single variable e, and parameters (α, ≳): $\mathit{F}\mathrm{(}\mathit{e,\alpha ,}\gtrsim \mathrm{)}\mathrm{=}\mathit{\alpha e}\mathrm{+}\frac{\gtrsim \mathit{e}}{\sqrt{\mathit{\alpha }\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\mathrm{)}\mathrm{+}{\gtrsim }^{\mathrm{2}}{\mathit{e}}^{\mathrm{2}}}}\mathrm{-}\mathrm{1}\mathrm{+}\mathit{\alpha }\mathrm{=}\mathrm{0.}$(35)Here F(e,a, ≳) is properly defined for (e,a, ≳) in the domain D_{e,α, ≳} defined by e ∈ [ 0,1 ] ,α ∈ ] 0,1 ] , ≳ ∈ ] 0, + ∞ [, as in this domain, 1 − e² + ≳ ²e²/α> 0. We also have $\frac{\mathit{\partial F}}{\mathit{\partial e}}\mathrm{(}\mathit{e,\alpha ,}\gtrsim \mathrm{)}\mathrm{=}\mathit{\alpha }\mathrm{+}\frac{\mathit{\alpha }\gtrsim }{{\left(\mathit{\alpha }\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\mathrm{)}\mathrm{+}{\gtrsim }^{\mathrm{2}}{\mathit{e}}^{\mathrm{2}}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}}\mathrm{·}$(36)Thus, $\frac{\mathit{\partial F}}{\mathit{\partial e}}\mathrm{(}\mathit{e,\alpha ,}\gtrsim \mathrm{)}\mathit{>}\mathrm{0}$ on the domain D_{e,α, ≳} and F(e,α, ≳) is strictly increasing with respect to e for e ∈ [ 0,1 ]. Moreover, as 0 <α< 1, $\mathit{F}\mathrm{(}\mathrm{0}\mathit{,\alpha ,}\gtrsim \mathrm{)}\mathrm{=}\mathrm{-}\mathrm{1}\mathrm{+}\mathit{\alpha }\mathit{<}\mathrm{0}\mathrm{;} \mathit{F}\mathrm{(}\mathrm{1}\mathit{,\alpha ,}\gtrsim \mathrm{)}\mathrm{=}\mathrm{2}\mathit{\alpha }\mathit{>}\mathrm{0}\mathrm{;}$(37)and $\mathit{F}\mathrm{(}{\mathit{e}}_{\mathrm{0}}\mathit{,\alpha ,}\gtrsim \mathrm{)}\mathrm{=}\frac{\gtrsim {\mathit{e}}_{\mathrm{0}}}{\sqrt{\mathit{\alpha }\mathrm{(}\mathrm{1}\mathrm{-}{{\mathit{e}}_{\mathrm{0}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{+}{\gtrsim }^{\mathrm{2}}{\mathit{e}}_{\mathrm{0}}^{\mathrm{2}}}}\mathit{>}\mathrm{0.}$(38)The equation of collision (35) thus always has a single solution e_c in the interval ] 0,min(1,e₀) [. This ensures that this critical value of e will fulfil the condition (30). The corresponding value of the relative AMD ${\mathit{C}}_{\mathit{c}}\mathrm{(}\mathit{\alpha ,}\gtrsim \mathrm{)}\mathrm{=}C\mathrm{\left(}{\mathit{e}}_{\mathit{c}}\mathit{,}{\mathit{e}}_{\mathit{c}}^{\mathrm{\prime }}\mathrm{\right)}$ is then obtained through (29).

3.2. Critical AMD C${{c}}_{}\mathrm{(}\mathit{\alpha ,}\gtrsim \mathrm{)}$

We have thus demonstrated that for a given pair of ratios of semi-major axes and masses, (α, ≳), there is always a unique critical value C_c(α, ≳) of the relative AMD C= C/ Λ′ which defines the AMD-stability. The system of two planets is AMD-stable if and only if $C\mathrm{=}\frac{\mathit{C}}{{\mathrm{\Lambda }}^{\mathrm{\prime }}}\mathit{<}{\mathit{C}}_{\mathit{c}}\mathrm{(}\mathit{\alpha ,}\gtrsim \mathrm{)}\mathit{.}$(39)The value of the critical AMD C_c(α, ≳) is obtained by computing first the critical eccentricity e_c(α, ≳) which is the unique solution of the collisional Eq. (35) in the interval [ 0,1 ]. The critical AMD is then ${\mathit{C}}_{\mathit{c}}\mathrm{(}\mathit{\alpha ,}\gtrsim \mathrm{)}\mathrm{=}C\mathrm{\left(}{\mathit{e}}_{\mathit{c}}\mathit{,}{\mathit{e}}_{\mathit{c}}^{\mathrm{\prime }}\mathrm{\right)}$ (Eq. (29)), where the critical value ${\mathit{e}}_{\mathit{c}}^{\mathrm{\prime }}$ is obtained from e_c through Eq. (28). It is important to note that the critical AMD, and thus the AMD-stability condition, depends only on (α, ≳).

3.3. Behaviour of the critical AMD

We now analyse the general properties of the critical AMD function C_c(α, ≳). As $\frac{\mathit{\partial F}}{\mathit{\partial e}}\mathrm{(}\mathit{e,\alpha ,}\gtrsim \mathrm{)}\mathit{>}\mathrm{0}$, we can apply the implicit function theorem to the domain D_{e,α, ≳}, which then ensures that the solution of the collision Eq. (35), e_c(α, ≳), is a continuous function of ≳ (and even analytic for ≳ ∈ ] 0, + ∞ [). Moreover, on D_{e,α, ≳},$\frac{\mathit{\partial F}}{\mathit{\partial }\gtrsim }\mathrm{(}\mathit{e,\alpha ,}\gtrsim \mathrm{)}\mathrm{=}\frac{\mathit{\alpha e}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\mathrm{)}}{{\left(\mathit{\alpha }\left(\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\right)\mathrm{+}{\gtrsim }^{\mathrm{2}}{\mathit{e}}^{\mathrm{2}}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}}\mathrm{\ge }\mathrm{0.}$(40)We also have$\frac{\mathit{\partial }{\mathit{e}}_{\mathit{c}}}{\mathit{\partial }\gtrsim }\mathrm{(}\mathit{\alpha ,}\gtrsim \mathrm{)}\mathrm{=}\mathrm{-}\frac{{\mathit{e}}_{\mathit{c}}\mathrm{(}\mathrm{1}\mathrm{-}{{\mathit{e}}_{\mathit{c}}^{\mathrm{2}}}^{\mathrm{)}}}{{\left(\mathit{\alpha }\left(\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathit{c}}^{\mathrm{2}}\right)\mathrm{+}{\gtrsim }^{\mathrm{2}}{\mathit{e}}_{\mathit{c}}^{\mathrm{2}}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{+}\gtrsim }\mathrm{\le }\mathrm{0}\mathit{,}$(41)and e_c(α, ≳) is a decreasing function of ≳. For any given values of the semi-major axes ratio α, and masses ≳, we can thus find the critical value C_c(α, ≳) which allows for a collision (39). For the critical value C_c(α, ≳), a single solution corresponds to the tangency condition (Fig. 1), and this solution is obtained at the critical value e_c(α, ≳) for the eccentricity of the orbit . The values of the critical relative AMD C_c(α, ≳) are plotted in Fig. 2 versus α, for different values of ≳. Deriving Eq. (29)with respect to ≳, one obtains $\begin{array}{ccc}\frac{\mathit{\partial }C}{\mathit{\partial }\gtrsim }& \mathrm{=}& \sqrt{\mathit{\alpha }}\left(\mathrm{1}\mathrm{-}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}}\right)\\ & & \mathrm{+}\hspace{0.17em}\gtrsim \sqrt{\mathit{\alpha }}\frac{\mathit{e}}{\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}}}\frac{\mathit{\partial e}}{\mathit{\partial }\gtrsim }\mathrm{+}\frac{{\mathit{e}}^{\mathrm{\prime }}}{\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{\prime }\mathrm{2}}}}\frac{\mathit{\partial }{\mathit{e}}^{\mathrm{\prime }}}{\mathit{\partial }\gtrsim }\mathrm{·}\end{array}$(42)Using the two relations (Eqs. (28), (34)), this reduces to $\frac{\mathit{\partial }{\mathit{C}}_{\mathit{c}}\mathrm{(}\mathit{\alpha ,}\gtrsim \mathrm{)}}{\mathit{\partial }\gtrsim }\mathrm{=}\sqrt{\mathit{\alpha }}\left(\mathrm{1}\mathrm{-}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathit{c}}^{\mathrm{2}}}\right)\mathit{>}\mathrm{0.}$(43)Thus C_c(α, ≳) strictly increases with ≳. In the same way $\frac{\mathit{\partial }{\mathit{C}}_{\mathit{c}}\mathrm{(}\mathit{\alpha ,}\gtrsim \mathrm{)}}{\mathit{\partial \alpha }}\mathrm{=}\frac{\gtrsim }{\mathrm{2}\sqrt{\mathit{\alpha }}}\left(\mathrm{1}\mathrm{-}\frac{\mathrm{(}\mathrm{1}\mathrm{+}{{\mathit{e}}_{\mathit{c}}^{\mathrm{2}}}^{\mathrm{)}}}{\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathit{c}}^{\mathrm{2}}}}\right)\mathit{<}\mathrm{0}$(44)and C_c(α, ≳) decreases with α.

Fig. 2 Values of Cc(α, ≳) vs. α for the different ≳ values for which an explicit expression of Cc(α, ≳) is obtained.

Now that the general behaviour of the critical AMD C_c(α, ≳) is known, we can specify its explicit expression in a few special cases that are displayed in Table 1. The computations and the higher order developments can be found in Appendix C.

A development of C_c(α, ≳) for α − → 1 can also be obtained (see Appendix C). With η = 1 − α, we have ${\mathit{C}}_{\mathit{c}}\mathrm{(}\mathit{\alpha ,}\gtrsim \mathrm{)}\mathrm{~}\frac{\mathit{\gamma }}{\mathit{\gamma }\mathrm{+}\mathrm{1}}\frac{{\mathit{\eta }}^{\mathrm{2}}}{\mathrm{2}}\mathrm{+}\mathrm{O}\mathrm{\left(}{\mathit{\eta }}^{\mathrm{3}}\mathrm{\right)}\mathit{.}$(45)We now present two applications of the AMD-stability, the formation of planetary systems, and a classification of observed planetary systems.

4. AMD model of planet formation

Once the disc has been depleted, the last phase of terrestrial planet formation begins with a disc composed of planetary embryos and planetesimals (Morbidelli et al. 2012). In numerical simulations of this phase, the AMD has been used as a statistical measure for comparison with the inner solar system (Chambers 2001). In this part, we recall the very simple model of embryo accretion described in Laskar (2000) interpreting the N-body dynamics as AMD-exchange.

4.1. Collisional evolution

So far we have not made special simplifications, and the model simply preserves the mass and the momentum in the barycentric reference frame. We make an additional assumption here: between collisions, the evolution of the orbits is similar to the evolution they would have in the averaged system in the presence of chaotic behaviour. The orbital parameters will thus evolve in a limited manner, with a stochastic diffusion process, bounded by the conservation of the total AMD. As shown in Sect. 2.3, during a collision the AMD decreases (Eq. (24)). We assume the collisions to be totally inelastic with perfect merging.Indeed, Kokubo & Genda (2010) and Chambers (2013) have shown that the detailed mechanisms of collisions, such as the possibilities of hit-and-run or fragmentation of embryos, barely change the final architecture of simulated systems. The total AMD of the system will thus be constant between collisions, and will decrease during collisions. On the other hand, the AMD for a particle is of the order of $\frac{\mathrm{1}}{\mathrm{2}}\mathit{m}\sqrt{\mathit{a}}{\mathit{e}}^{\mathrm{2}}$. As the mass of the particle increases, its excursion in eccentricity will be more limited, and fewer collisions will be possible. The collisions will stop when the total AMD of the system is too small to allow for planetary collisions.

In the accretion process, we consider a planetesimal of semi-major axis a and its immediate neighbour, defined as the planetesimal with semi-major axis a′, such that there is no other planetesimal with semi-major axis between a and a′. In this case we can assume that α is close to 1 and, as explained in Appendix C, we use an approximation of the critical AMD value C_c(α,γ): ${\mathit{C}}_{\mathit{c}\mathrm{1}}\mathrm{(}\mathit{\alpha ,}\gtrsim \mathrm{)}\mathrm{=}\mathit{k}\mathrm{(}\gtrsim \mathrm{)}{\left(\frac{{\mathit{\delta }}_{\mathit{a}}}{\mathit{a}}\right)}^{\mathrm{2}}\mathit{,}$(46)where δ_a = a′ − a and k( ≳ ) = ≳ / (2( ≳ + 1)).

4.2. AMD-stable planetary distribution

In this section, we search for the laws followed by the planetary distribution of a model formed following the above assumptions. We thus start from an arbitrary distribution of mass of planetesimals ρ(a), and let the system evolve under the previous rules. We search for the condition of AMD-stable planetary systems, obtained by random accretion of planetesimals. This condition requires that the final AMD value cannot allow for planetary crossing among the planets. Let C be the value of the AMD at the end of the accretion process. If we consider the planetesimal of semi-major axis a, its mass will continue to increase by accretion with a body of semi-major axis a′>a, as long as $C\mathrm{=}\frac{\mathit{C}}{{\mathrm{\Lambda }}^{\mathrm{\prime }}}\mathrm{\ge }{\mathit{C}}_{\mathit{c}\mathrm{1}}\mathrm{=}\mathit{k}{\left(\frac{{\mathit{\delta }}_{\mathit{a}}}{\mathit{a}}\right)}^{\mathrm{2}}\mathrm{·}$(47)The initial linear density of mass is ρ(a). As a′ is the closest neighbour to a, we can assume that all the planetesimals initially between a and a′ have been absorbed by the two bodies of mass m(a) and m(a′). At first order with respect to δ_a/a, we have thus $\mathit{m}\mathrm{\left(}{\mathit{a}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{~}\mathit{m}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{~}\mathit{\rho }\mathrm{\left(}\mathit{a}\mathrm{\right)}{\mathit{\delta }}_{\mathit{a}}\mathit{,}$(48)and from (47) at the limiting case we have $\frac{\begin{array}{c}\mathrm{˜}\\ \mathit{C}\end{array}}{{\mathit{\delta }}_{\mathit{a}}\mathit{\rho }\mathrm{\left(}\mathit{a}\mathrm{\right)}\sqrt{\mathit{a}}}\mathrm{=}\mathit{k}{\left(\frac{{\mathit{\delta }}_{\mathit{a}}}{\mathit{a}}\right)}^{\mathrm{2}}\mathit{,}$(49)where $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{C}\end{array}\mathrm{=}\mathit{C}\mathit{/}\sqrt{\mathit{\mu }}$. We have thus ${\mathit{\delta }}_{\mathit{a}}\mathrm{=}{\left(\frac{\begin{array}{c}\mathrm{˜}\\ \mathit{C}\end{array}}{\mathit{k}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}{\mathit{a}}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{\rho }}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{3}}\mathit{,}$(50)and from Eq. (48) $\mathit{m}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{=}{\left(\frac{\begin{array}{c}\mathrm{˜}\\ \mathit{C}\end{array}}{\mathit{k}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}{\mathit{a}}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{\rho }}^{\mathrm{2}\mathit{/}\mathrm{3}}\mathit{.}$(51)

4.3. Scaling laws with initial mass distribution ρ(a) ∝ a^p

Using the previous relations, we can compute the resulting systems for various initial mass distribution, in particular for ρ(a) = ζa^p. From Eq. (51), we obtain for two consecutive planets $\frac{\mathit{m}}{{\mathit{m}}^{\mathrm{\prime }}}\mathrm{=}{\mathit{\alpha }}^{\mathrm{1}\mathit{/}\mathrm{2}}{\left(\frac{\mathit{\rho }\mathrm{\left(}\mathit{a}\mathrm{\right)}}{\mathit{\rho }\mathrm{\left(}{\mathit{a}}^{\mathrm{\prime }}\mathrm{\right)}}\right)}^{\mathrm{2}\mathit{/}\mathrm{3}}\mathrm{=}{\mathit{\alpha }}^{\frac{\mathrm{1}}{\mathrm{2}}\mathrm{+}\frac{\mathrm{2}}{\mathrm{3}}\mathit{p}}$(52)and from Eq. (46), as lim_{α − → 1} ≳ = 1, we thus have $\mathit{k}\mathrm{(}\gtrsim \mathrm{)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}}$ in all cases. Relation (50) can be rewritten ${\mathit{\delta }}_{\mathit{a}}{\mathit{a}}^{\frac{\mathit{p}}{\mathrm{3}}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{=}\mathrm{(}\mathrm{4}\begin{array}{c}\mathrm{˜}\\ \mathit{C}\end{array}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{3}}{\mathit{\zeta }}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{3}}\mathit{\delta n,}$(53)where δn = 1 is the increment from planet a to a′. By integration, this difference relation becomes for $\mathit{p}\ne \mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}$, ${\mathit{a}}^{\frac{\mathrm{2}\mathit{p}\mathrm{+}\mathrm{3}}{\mathrm{6}}}\mathrm{=}{\mathit{a}}_{\mathrm{0}}^{\frac{\mathrm{2}\mathit{p}\mathrm{+}\mathrm{3}}{\mathrm{6}}}\mathrm{+}\frac{\mathrm{2}\mathit{p}\mathrm{+}\mathrm{3}}{\mathrm{6}}{\left(\frac{\mathrm{4}\begin{array}{c}\mathrm{˜}\\ \mathit{C}\end{array}}{\mathit{\zeta }}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathit{n}\mathit{.}$(54)For $\mathit{p}\mathrm{=}\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}$, we obtain $\log \mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{=}\log \mathrm{\left(}{\mathit{a}}_{\mathrm{0}}\mathrm{\right)}\mathrm{+}\mathrm{(}\mathrm{4}\begin{array}{c}\mathrm{˜}\\ \mathit{C}\end{array}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{3}}{\mathit{\zeta }}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{3}}\mathit{n}\mathit{.}$(55)In particular, for p = 0 (constant distribution and ρ(a) = ζ), from Eq. (53), we have $\sqrt{\mathit{a}}\mathrm{=}\sqrt{{\mathit{a}}_{\mathrm{0}}}\mathrm{+}{\left(\frac{\begin{array}{c}\mathrm{˜}\\ \mathit{C}\end{array}}{\mathrm{2}\mathit{\zeta }}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathit{n}\mathit{.}$(56)For the masses, from Eq. (51), we have for large n (or if a₀ is small) $\mathit{m}\mathrm{\left(}\mathit{n}\mathrm{\right)}\mathrm{~}\mathrm{(}\mathrm{2}{\begin{array}{c}\mathrm{˜}\\ \mathit{C}\end{array}}^{\mathrm{2}}\mathit{\zeta }{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathit{n}\mathit{.}$(57)For p = − 3 / 2, we find a power law similar to the Titius-Bode law for planetary distribution². The different expressions deduced from this model of planetary accretion are listed in Table 2.

The previous analytical results were tested on a numerical model of our accretion scheme (Laskar 2000). The model was designed to fulfil the conditions (Eqs. (12), (13)) of Sect. 2.3. Five thousand simulations were started with a large number of orbits (10 000) and followed in order to look for orbit intersections. When an intersection occurs, the two bodies merge into a new one whose orbital parameters are determined by the collisional Eqs. (12), (13) (see Appendix A). Between collisions, the orbits do not evolve, apart from a diffusion of their eccentricities, which fulfils the condition of conservation of the total AMD. This is roughly what would occur in a chaotic secular motion.

The main parameter of these simulations is the final AMD value, C. Because the AMD decreases during collisions, and in order to obtain final systems with a given value C of the AMD, the eccentricities were increased by a small amount in order to raise the AMD to the desired final value. This is justified as close encounters can also increase the AMD value. Indeed, N-body simulations (Chambers 2001; Raymond et al. 2006) present a first phase of AMD increase at the beginning of the simulations before the AMD decreases and converges to its final value. These simulations were extremely rapid to integrate as the orbital motion of the orbits was not integrated. Instead, we looked here for collisions of ellipses which fulfil the conservation of mass and of linear momentum. These simulations were thus started with a large number of initial bodies (10 000) and continued until their final evolution. The different runs resulted in different numbers of planets, which ranged from four to nine, but the averaged values fitted very well with the analytical results of Table 2 (Laskar 2000).

5. AMD-classification of planetary systems

Here we show how the AMD-stability can be used as a criterion to derive a classification of planetary systems. In Sect. 3.2 we saw that in the secular approximation, the stability of a pair of planets is determined by the computation of $\mathit{\beta }\mathrm{=}\frac{C}{{\mathit{C}}_{\mathit{c}}}\mathrm{=}\frac{\mathit{C}}{{\mathrm{\Lambda }}^{\mathrm{\prime }}{\mathit{C}}_{\mathit{c}}}\mathrm{·}$(58)We call β the AMD-stability coefficient. For pairs of planets, β< 1 means that collisions are not possible. The pairs of planets are then called AMD-stable. We naturally extend this definition to multiple planet systems. A system is AMD-stable if every adjacent pair is AMD-stable³. We can also define an AMD-stability coefficient regarding the collision with the star. We define β_S, the AMD-stability coefficient of the pair formed by the star and the innermost planet. For this pair, we have α = 0 and thus C_c = 1. With this simplification β_S = C/ Λ, where Λ is the circular momentum of the innermost planet.

5.1. Sample studied and methods of computation

We have studied the AMD-stability of some systems referenced in the exoplanet.eu catalogue (Schneider et al. 2011). From the catalogue, we selected 131 systems that have measured masses, semi-major axes, and eccentricities for all their planets. Since the number of systems with known mutual inclinations is too small, we assumed the systems to be almost coplanar. This claim is supported by previous statistical studies that constrain the observed mutual inclinations distribution (Fang & Margot 2012; Lissauer et al. 2011; Fabrycky et al. 2014; and Figueira et al. 2012). For some systems where the uncertainties were not provided, we consulted the original papers or the Exoplanet Orbit Database⁴ (Wright et al. 2011).

Fig. 3 Cumulative distributions of the period ratios of adjacent planets in the sample studied here, of the AMD-stable systems (both weak and strong), AMD-unstable systems, and for all the systems referenced in the exoplanet.eu database.

We compare the cumulative distribution of the adjacent planets’ period ratios in our sample and that of all the multiplanetary systems in the database exoplanet.eu (see Fig. 3). The sample is biased toward higher period ratios. Indeed, most of the multiplanetary systems in the database come from the Kepler data. Since these systems are mostly tightly packed ones, their period ratios are rather small. However, the majority of them do not have measured eccentricities and are consequently excluded from this study. Our sample thus mostly contains systems detected by radial velocities (RV) methods that have, on average, higher period ratios.

Since all the AMD computations are done with the relative quantities α and γ, we can use equivalent quantities that are measured more precisely in observations than the masses and semi-major axis. We used the period ratios elevated to the power 3/2 instead of the semi-major axis ratios, and the minimum mass msin(i) for RV system. This is not a problem for the computation of γ; if we assume that the systems are close to coplanarity, then $\mathit{\gamma }\mathrm{=}\frac{\mathit{m}}{{\mathit{m}}^{\mathrm{\prime }}}\mathrm{=}\frac{\mathit{m}\sin \mathrm{\left(}\mathit{i}\mathrm{\right)}}{{\mathit{m}}^{\mathrm{\prime }}\sin \mathrm{\left(}\mathit{i}\mathrm{\right)}}\mathrm{\simeq }\frac{\mathit{m}\sin \mathrm{\left(}\mathit{i}\mathrm{\right)}}{{\mathit{m}}^{\mathrm{\prime }}\sin \mathrm{\left(}{\mathit{i}}^{\mathrm{\prime }}\mathrm{\right)}}\mathrm{·}$(59)Even though we assume the systems to be coplanar, we want to take into account the contribution of mutual inclinations to the AMD. Since we only have access to the eccentricities, we define the coplanar AMD of a system C_p, as the AMD of the same system if it was coplanar. We can also define ${\mathit{\beta }}_{\mathit{p}}\mathrm{=}{\mathit{C}}_{\mathit{p}}\mathit{/}\mathrm{\left(}{\mathrm{\Lambda }}^{\mathrm{\prime }}{\mathit{C}}_{\mathit{c}}\mathrm{\right)}\mathit{,}$(60)which is the coplanar AMD-stability coefficient. Motivated by both theoretical arguments on chaotic diffusion in the secular dynamics (Laskar 1994, 2008) and observed correlations in statistical distributions (Xie et al. 2016), we assume that the AMD contribution of mutual inclinations is equal to that of the eccentricities. This hypothesis is equivalent to assuming the average equipartition of the AMD among the secular degrees of freedom for a chaotic system. The total AMD is thus twice the measured AMD, and in this study we use $\mathit{\beta }\mathrm{=}\mathrm{2}{\mathit{\beta }}_{\mathit{p}}\mathit{.}$(61)We can also define a coplanar AMD-stability coefficient associated with the star, and similarly we set β_S=2β_Sp. We then compute the coefficients β_S and β for each pair and the associated uncertainty distributions. We list the results of the analysis in Table 3. In the considered dataset, 70 systems are AMD-stable. The majority of the highest multiplicity systems are AMD-unstable. In Fig. 4 we plot the probability distribution of β for the considered systems.

Fig. 4 Probability distribution function of β for the sample studied. The systems are grouped by multiplicity. The vertical line β = 1 marks the separation between AMD-stable and AMD-unstable pairs.

5.2. Propagation of uncertainties

The uncertainties are propagated using Monte Carlo simulations of the distributions. After determining the distributions from the input quantities (m,a,e), we generate 10 000 values for each of these parameters. We then compute the derived quantities (α,γ,C_c,β...) in these 10 000 cases and the associated distributions.

For masses (or msin(i) if no masses are provided) and periods, we assume a Gaussian uncertainty centred on the value referenced in the database and with standard deviation, the half width of the confidence interval. The distributions are truncated to 0.

For eccentricity distributions, the previous method does not provide satisfying results. Most of the Gaussian distributions constructed with the mean value and confidence interval given in the catalogue make negative eccentricities probable (in the case of almost circular planets with a large upper bound on the eccentricity). One solution is to assume that the rectangular eccentricity coordinates (ecosω,esinω) are Gaussian. Since the average value of ω has no importance in the computation of the eccentricity distribution, we assume it to be 0. Therefore, esinω has 0 mean. We define the distribution of $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{e}\end{array}\mathrm{=}\mathit{e}\cos \mathit{\omega }$ as a Gaussian distribution with the mean value, the value referenced in the catalogue, and standard deviation, the half-width of the confidence interval. If we assume esinω has the same standard deviation as ecosω, we have $\mathit{e}\sin \mathit{\omega }\mathrm{=}\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{e}\end{array}\mathrm{-}\mathrm{⟨}\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{e}\end{array}\mathrm{⟩}$. The distribution of e is then deduced from that of $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{e}\end{array}$ using $\mathit{e}\mathrm{=}\sqrt{{\begin{array}{c}\mathrm{˜}\\ \mathit{e}\end{array}}^{\mathrm{2}}\mathrm{+}\mathrm{(}\begin{array}{c}\mathrm{˜}\\ \mathit{e}\end{array}\mathrm{-}\mathrm{⟨}\begin{array}{c}\mathrm{˜}\\ \mathit{e}\end{array}\mathrm{⟩}{\mathrm{)}}^{\mathrm{2}}}\mathit{.}$(62)Using the Gaussian assumption means that some masses or periods can take values close to 0 with a small probability. This causes the distributions of α or γ to diverge if it happens that a′ or m′ can take values close to 0. To address this issue, a linear expansion around the mean value is used for the quotients, for example for α, $\mathit{\alpha }\mathrm{=}\frac{\mathit{a}}{{\mathit{a}}^{\mathrm{\prime }}}\mathrm{=}\frac{\mathit{a}}{\mathrm{⟨}{\mathit{a}}^{\mathrm{\prime }}\mathrm{⟩}}\left(\mathrm{1}\mathrm{-}\frac{\mathrm{\Delta }{\mathit{a}}^{\mathrm{\prime }}}{\mathrm{⟨}{\mathit{a}}^{\mathrm{\prime }}\mathrm{⟩}}\right)\mathit{,}$(63)with Δa′ = a′ − ⟨ a′ ⟩.

Fig. 5 Architecture of the strong AMD-stable systems. Each planet is represented by a circle whose size is proportional to log 10(m), where m is the mass of the planet. The colour represents the AMD-stability coefficient of the inner pair associated with the planet (in particular, the first planet is represented with the AMD-stability coefficient associated with the star βS). This means that a red planet can collide with its inner neighbour.

Fig. 6 Architecture of the weak AMD-stable systems. Each planet is represented by a circle whose size is proportional to log 10(m), where m is the mass of the planet. The colour represents the AMD-stability coefficient of the inner pair associated with the planet (in particular, the first planet is represented with the AMD-stability coefficient associated with the star βS).

5.3. AMD-stable systems

As said above, we consider AMD-stable a system where collisions between planets are impossible because of the dynamics ruled by AMD (maxβ< 1). In addition, if the AMD-stability coefficient of the star β_S< 1 (resp. β_S> 1), the system is defined as strong AMD-stable (resp. weak AMD-stable).

5.3.1. Strong AMD-stable systems

The strong AMD-stable systems can be considered dynamically stable in the long term. In Fig. 5, we plot the architecture of the strong AMD-stable systems. If the system is out of the mean motion resonance (MMR) islands, the AMD will not increase and therefore no collision between planets or with the star can occur. We can see in Fig. 3 that the AMD-stable systems have period ratios on average larger than those from the sample.

5.3.2. Weak AMD-stable systems

As defined above, in a weak AMD-stable system, no planetary collisions can occur, but the innermost planet can increase its eccentricity up to 1 and collide with the star. We separate these systems from the strong AMD-stable ones because the system can still lose a planet only by AMD exchange. However, the remaining system will not be affected by the destruction of the inner planet. In Fig. 6, we plot their architecture. In these systems, the inner planet is much closer to the star than the others.

5.4. AMD-unstable systems

The AMD-unstable systems have at least one unstable planet pair, but as we show in Fig. 7 where we plot the architecture of these systems, this category is not homogeneous. It gathers high multiplicity systems where planets are too close to each other for the criterion to be valid, multiscale systems like the solar system where the inner system is AMD-unstable owing to its small mass in comparison to the outer part, systems experiencing mean motion resonances, etc. In all these cases, an in-depth dynamical study is necessary to determine the short- or long-term stability of these systems. In the following sections, we detail how the AMD-stability and AMD driven dynamics can help to understand these systems.

5.4.1. Hierarchically AMD-stable systems (solar system-like)

We first consider the solar system. Owing to the large amount of AMD created by the giant planets, the inner system is not AMD-stable (Laskar 1997). However, long-term secular and direct integrations have shown that the inner system has a probability of becoming unstable of only 1% over 5 Gyr (Laskar 2008; Laskar & Gastineau 2009; Batygin et al. 2015). Indeed, the chaotic dynamics is mainly restricted to the inner system, while the outer system is mostly quasi-periodic. However, when AMD exchange occurs between the outer and inner systems, the inner system becomes highly unstable and can lose one or several planets.

In Fig. 8 we plot the AMD and the circular momenta of the solar system planets. We see that the inner system planets (resp. the outer planets) have comparable AMD values. Laskar (2008) showed the absence of exchange between the two parts of the system. Moreover, the spacing of the planets follows surprisingly well the distribution laws mentioned in Sect. 4 if we consider the two parts of the system separately (Laskar 2000).

We see that given an analysis of the secular dynamics, the tools developed above can still help to understand how an a priori unstable system can still exist in its current state. The case of AMD-unstable systems with an AMD-stable part is not restricted to the solar system. If we look for systems where the outer part is AMD-stable while the whole system is AMD-unstable, we find four similar systems in our sample as shown in Fig. 9. We call these systems hierarchic AMD-stable systems.

5.4.2. Resonant systems

As shown by the cumulative distributions plotted in Fig. 3, the AMD-unstable systems have period ratios that are lower than the AMD-stable systems. For example, in our sample 66% of the adjacent pairs of the AMD-unstable systems have a period ratio below 4, whereas this proportion is only 33% for the AMD-stable systems. For period ratios close to small integer ratios, the MMR can rule a great part of the dynamics. Particularly, if pairs of planets are trapped in a large resonance island, the system could be dynamically stable even if it is not AMD-stable.

Individual dynamical studies are necessary in order to claim that a system is stable due to a MMR. We note, however, that the AMD-unstable systems considered here are statistically constituted of more resonant pairs than a typical sample in the catalogue.

5.4.3. Concentration around MMR

Here we want to test whether the unstable systems have been drawn randomly from the exoplanet.eu catalogue or if the period ratios of the pairs of adjacent planets are statistically closer to small integer ratios. We denote ℋ₀, the hypothesis that the period ratios of the unstable systems have been drawn randomly from the catalogue distribution. We consider only the period ratios lower than 4 because the higher ones are not significant for a study of the MMR. We plot in Fig. 10 the cumulative distribution of the period ratios of the AMD-unstable systems, as well as the cumulative distribution of the period ratios of all the systems of exoplanet.eu. We call R_u the set of period ratios of the AMD-unstable planets.

Fig. 7 Architecture of the AMD-unstable systems. Each planet is represented by a circle whose size is proportional to log 10(m), where m is the mass of the planet. The colour represents the AMD-stability coefficient of the inner pair associated with the planet (in particular, the first planet is represented with the AMD-stability coefficient associated with the star βS).

We first test ℋ₀ via a Kolmogorov-Smirnov test (Lehmann & Romano 2006) between the sample R_u and the period ratios of the catalogue. The test fails to reject ℋ₀ with a p-value of about 9%. Therefore, we cannot reject ℋ₀ on the global shape of the distribution of the AMD-unstable period ratios.

However, we want to determine if the AMD-unstable pairs are close to small integer ratios, which means studying the fine structure of the distribution. We propose here another method for demonstrating this.

Let us denote the probability density of the catalogue period ratios f and the associated cumulative distribution ℱ. Let us consider an interval I_x = (x,x + Δx), if we assume ℋ₀, the probability for a ratio r to be in I_x is $\begin{array}{ccc}\mathrm{𝒫}\mathrm{(}\mathit{r}\mathrm{\in }{\mathit{I}}_{\mathit{x}}\mathrm{|}{\mathrm{\mathscr{H}}}_{\mathrm{0}}\mathrm{)}& \mathrm{=}& \int \begin{array}{c}\mathit{x}\mathrm{+}\mathrm{\Delta }\mathit{x}\\ \mathit{x}\end{array}\mathit{f}\mathrm{\left(}{\mathit{r}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{d}{\mathit{r}}^{\mathrm{\prime }}\\ & \mathrm{=}& \mathrm{ℱ}\mathrm{(}\mathit{x}\mathrm{+}\mathrm{\Delta }\mathit{x}\mathrm{)}\mathrm{-}\mathrm{ℱ}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\mathrm{\Delta ℱ}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathit{.}\end{array}$(64)Therefore, the probability that more than k₀ pairs out of N pairs drawn randomly from the distribution ℱ are in I_x is $\mathrm{𝒫}\mathrm{\left(}\mathit{x,}{\mathit{k}}_{\mathrm{0}}\mathrm{\right|}{\mathrm{\mathscr{H}}}_{\mathrm{0}}\mathrm{)}\mathrm{=}\sum _{\mathit{k}\mathrm{=}{\mathit{k}}_{\mathrm{0}}}^{\mathit{N}}\left(\begin{array}{c}\mathit{N}\\ \mathit{k}\end{array}\right)\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{\Delta ℱ}{\mathrm{)}}^{\mathit{N}\mathrm{-}\mathit{k}}\mathrm{\Delta }{\mathrm{ℱ}}^{\mathit{k}}\mathit{.}$(65)This is the probability of a binomial trial. From now on, N designates the number of period ratios in R_u.

For a given Δx, we can compute the probability , where k_u(x) is the number of pairs from R_u in I_x. This probability measures the likelihood of a concentration of elements of R_u around x, assuming they were drawn from ℱ. We choose Δx = 0.05 and plot the function P(x,R_u) as well as k_u in Fig. 11. We observe that the concentrations around 3 and 2 are very unlikely with a probability of 2.3 × 10^-3 and 2.9 × 10^-3, respectively. Other peaks appear for x = 1.4 and around 3.8. However, P(x,R_u) is the probability of seeing a concentration around x precisely. It is not the probability of observing a concentration somewhere between 1 and 4.

To demonstrate that the concentrations around 2 and 3 are meaningful, we compare this result to samples drawn randomly from the distribution ℱ. We create 10,000 datasets R_α by drawing N points from ℱ and compute P(x,R_α) for these datasets. Then, we record the two minimum values (at least distant by more than Δx) and plot the distribution of these minima in Fig. 12. From these simulations, we see that on average 17.2% of the samples have a minimum as small as R_u. However, the presence of a secondary peak as significant as the second one of R_u has a probability of 1.3%. Moreover, the R_u concentrations are clearly situated around low integer ratios, which would not be the case in general for a randomly generated sample.

We thus demonstrated here that the AMD-unstable systems period ratios are significantly concentrated around low integer ratios. While we do not prove that the pairs of planets close to these ratios are actually in MMR, this result further motivates investigations toward the behaviour of these pairs in a context of secular chaotic dynamics.

Fig. 8 Circular momenta and AMD of the solar system planets. The horizontal lines correspond to the respective mean Λk and Ck for the inner and outer system.

Fig. 9 Examples of systems where the outer part is AMD-stable and the inner part becomes AMD-stable if considered alone.

Fig. 10 Cumulative period ratios of the AMD-unstable systems and of all systems in the catalogue. The vertical line represents low order MMR.

Fig. 11 Red curve: probability of observing more than k0 pairs in the interval Ix given a probability density f. Green curve: k0 is the number of pairs in Ix for the unstable systems Ru.

Fig. 12 Probability density function of the two first minimum values of P(x,Rα) from 10 000 samples Rα drawn from the ℱ distribution.

6. Conclusion

The angular momentum deficit (AMD; Laskar 1997, 2000) is a key parameter for the understanding of the outcome of the formation processes of planetary systems (e.g. Chambers 2001; Morbidelli et al. 2012; Tremaine 2015). We have shown here how AMD can be used to derive a well-defined stability criterion: the AMD-stability. The AMD-stability of a system can be checked by the computation of the critical AMD C_c (Eq. (35)) and AMD-stability coefficients β_i that depend only on the eccentricities and ratios of semi-major axes and masses (Eqs. (29), (35), (39), (58)). This criterion thus does not depend on the degeneracy of the masses coming from radial velocity measures. Moreover, the uncertainty on the relative inclinations can be compensated by assuming equipartition of the AMD between eccentricities and inclinations.

AMD-stability will ensure that in the absence of mean motion resonances, the system is long-term stable. A rapid estimate of the stability of a system can thus be obtain by a short-term integration and the simple computation of the AMD-stability coefficients.

We have also proposed here a classification of the planetary systems based on AMD-stability (Sect. 5). The strong AMD-stable systems are the systems where no planetary collisions are possible, and no collisions of the inner planet with its central star, while the weak AMD-stable systems allow for the collision of the inner planet with the central star. The AMD-unstable systems are the systems for which the AMD-coefficient does not prevent the possibility of collisions. The solar system is AMD-unstable, but it belongs to the subcategory of hierarchical AMD-stable systems that are AMD-unstable but become AMD-stable when they are split into two parts (giant planets and terrestrial planets for the solar system) (Laskar 2000). Out of the 131 studied systems from exoplanet.eu, we find 48 strong AMD-stable, 22 weak AMD-stable, and 61 AMD-unstable systems, including 5 hierarchical AMD-stable systems.

As for the solar system, the AMD-unstable systems are not necessarily unstable, but determining their stability requires some further dynamical analysis. Several mechanisms can stabilise AMD-unstable systems. The absence of secular chaotic interactions between parts of the systems, like in the solar system case, or the presence of mean motion resonances, protecting pairs of planets from collision. In this case, the AMD-stability classification is still useful in order to select the systems that require this more in-depth dynamical analysis. It should also be noted that the discovery of additional planets in a system will require a revision of the computation of the AMD-stability of the system. This additional planet will always increase the total AMD, and thus the maximum AMD-coefficient of the system, decreasing its AMD-stability unless it is split into two subsystems.

In the present work, we have not taken into account mean motion resonances (MMR) and the chaotic behaviour resulting from their overlap. We aim to take these MMR into consideration in a forthcoming work. Indeed, the criterion regarding MMR developed by Wisdom (1980) or more recently Deck et al. (2013) may help to improve our stability criterion by considering the MMR chaotic zone as a limit for stability instead of the limit considered here that is given by the collisions of the orbits (Eqs. (35)). The drawback will then most probably be giving up the rigorous results that we have established in Sect. 3, and allowing for additional empirical studies. The present work will in any case remain the backbone of these further studies. Note: The AMD-stability coefficients of selected planetary systems are available on the exoplanet.eu database and at the CDS.

References

Appendix A: Intersection of planar orbits

Fig. A.1 Elliptical orbit, as the intersection of the cone r2 = x2 + y2 and the plane P·r + r = p.

In this section, we present an efficient algorithm for the computation of the intersection of two elliptical orbits in the plane, following (Albouy 2002). Let us consider an elliptical orbit defined by (μ,r,ṙ) and let G = r ∧ ṙ be the angular momentum per unit of mass. The Laplace vector ${P}\mathrm{=}\frac{{ṙ}\mathrm{\wedge }{G}}{\mathit{\mu }}\mathrm{-}\frac{{r}}{\mathit{r}}$(A.1)is an integral of the motion with coordinates (ecosω,esinω). One has ${P}\mathrm{·}{r}\mathrm{=}\frac{{\mathit{G}}^{\mathrm{2}}}{\mathit{\mu }}\mathrm{-}\mathit{r}\mathrm{=}\mathit{p}\mathrm{-}\mathit{r,}$(A.2)where p = a(1 − e²) is the parameter of the ellipse. Let r = (x,y) in the plane. We can consider the ellipse in three-dimensional space (x,y,r) as the intersection of the cone ${\mathit{r}}^{\mathrm{2}}\mathrm{=}{\mathit{x}}^{\mathrm{2}}\mathrm{+}\mathit{y}2$(A.3)with the plane defined by Eq. (A.2), that is $\mathit{x}\mathrm{\left(}\mathit{e}\cos \mathit{\omega }\mathrm{\right)}\mathrm{+}\mathit{y}\mathrm{\left(}\mathit{e}\sin \mathit{\omega }\mathrm{\right)}\mathrm{+}\mathit{r}\mathrm{=}\mathit{p}\mathit{.}$(A.4)If we consider now two orbits and . Their intersection is easily obtained as the intersection of the line of intersection of the two planes $\begin{array}{ccc}& & \mathit{x}\mathrm{\left(}{\mathit{e}}_{\mathrm{1}}\cos {\mathit{\omega }}_{\mathrm{1}}\mathrm{\right)}\mathrm{+}\mathit{y}\mathrm{\left(}{\mathit{e}}_{\mathrm{1}}\sin {\mathit{\omega }}_{\mathrm{1}}\mathrm{\right)}\mathrm{+}\mathit{r}\mathrm{=}{\mathit{p}}_{\mathrm{1}}\\ & & \mathit{x}\mathrm{\left(}{\mathit{e}}_{\mathrm{2}}\cos {\mathit{\omega }}_{\mathrm{2}}\mathrm{\right)}\mathrm{+}\mathit{y}\mathrm{\left(}{\mathit{e}}_{\mathrm{2}}\sin {\mathit{\omega }}_{\mathrm{2}}\mathrm{\right)}\mathrm{+}\mathit{r}\mathrm{=}{\mathit{p}}_{\mathrm{2}}\end{array}$(A.5)with the cone of equation r² = x² + y². Depending on the initial conditions, if and are distinct, we will get either 0, 1, or 2 solutions.

Appendix B: AMD in the averaged equations

In this section we show that the AMD conserves the same form in the averaged planetary Hamiltonian at all orders. More generally, this is true for any integral of H which does not depend on the longitude λ_k. Let $\mathit{H}\mathrm{=}{\mathit{H}}_{\mathrm{0}}\mathrm{\left(}\mathrm{\Lambda }\mathrm{\right)}\mathrm{+}\mathit{\epsilon }{\mathit{H}}_{\mathrm{1}}\mathrm{\left(}\mathrm{\Lambda }\mathit{,\lambda ,J,\phi }\mathrm{\right)}$(B.1)be a perturbed Hamiltonian system, and let K(Λ,J,φ) be a first integral of H(Λ,λ,J,φ) (such that { K,H } = 0, where { ·,· } is the usual Poisson bracket), and independent of λ. In addition, let ${\mathit{H}}^{\mathrm{\prime }}\mathrm{=}{\mathrm{e}}^{{\mathit{L}}_{\mathit{S}}}\mathit{H}$(B.2)be a formal averaging of H with respect to λ. If S(Λ,λ,J,φ) is a generator defined as below (Eqs. (B.4), (B.6), (B.7)), such that H′ is independent of λ. Then, K is an integral of H′, i.e. $\mathrm{\left\{}\mathit{K,}{\mathit{H}}^{\mathrm{\prime }}\mathrm{\right\}}\mathrm{=}\mathrm{0.}$(B.3)The generator S = εS₁ + ε²S₂ + ε³S₃ +··· is obtained formally through the identification order by order $\begin{array}{ccc}& & {\mathit{H}}_{\mathrm{0}}^{\mathrm{\prime }}\mathrm{=}{\mathit{H}}_{\mathrm{0}}\\ & & {\mathit{H}}_{\mathrm{1}}^{\mathrm{\prime }}\mathrm{=}{\mathit{H}}_{\mathrm{1}}\mathrm{+}\mathrm{\left\{}{\mathit{H}}_{\mathrm{0}}\mathit{,}{\mathit{S}}_{\mathrm{1}}\mathrm{\right\}}\\ & & {\mathit{H}}_{\mathrm{2}}^{\mathrm{\prime }}\mathrm{=}\mathrm{\left\{}{\mathit{H}}_{\mathrm{0}}\mathit{,}{\mathit{S}}_{\mathrm{2}}\mathrm{\right\}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\left\{}\mathrm{\right\{}{\mathit{H}}_{\mathrm{0}}\mathit{,}{\mathit{S}}_{\mathrm{1}}\mathrm{\left\}}\mathit{,}{\mathit{S}}_{\mathrm{1}}\mathrm{\right\}}\mathrm{+}\mathrm{\left\{}{\mathit{H}}_{\mathrm{1}}\mathit{,}{\mathit{S}}_{\mathrm{1}}\mathrm{\right\}}\\ & & \end{array}$(B.4)For any function G(Λ,λ,J,φ), let $\mathrm{⟨}\mathit{G}{\mathrm{⟩}}_{{\lambda }}\mathrm{=}\frac{\mathrm{1}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathit{m}}}\mathrm{\int }\mathit{G}{\mathrm{d}}^{\mathit{m}}{\lambda }$(B.5)be the average of G over all the angles λ_k. For each n ≥ 1, the S_n is obtained through the resolution of an equation of the form ${\mathit{H}}_{\mathit{n}}^{\mathrm{\prime }}\mathrm{=}\mathrm{\left\{}{\mathit{H}}_{\mathrm{0}}\mathit{,}{\mathit{S}}_{\mathit{n}}\mathrm{\right\}}\mathrm{+}{\mathrm{ℛ}}_{\mathit{n}}\mathit{,}$(B.6)where ℛ_n belongs to ℒ(H₀,H₁,S₁,...,S_{n − 1}), the Lie algebra generated by (H₀,H₁,S₁,...,S_{n − 1}). ${\mathit{H}}_{\mathit{n}}^{\mathrm{\prime }}$ will be the averaged part of ℛ_n, ⟨ ℛ_n ⟩ _λ, and S_n is obtained by solving the homological equation $\mathrm{\left\{}{\mathit{H}}_{\mathrm{0}}\mathit{,}{\mathit{S}}_{\mathit{n}}\mathrm{\right\}}\mathrm{=}\mathrm{⟨}{\mathrm{ℛ}}_{\mathit{n}}{\mathrm{⟩}}_{{\lambda }}\mathrm{-}{\mathrm{ℛ}}_{\mathit{n}}\mathit{.}$(B.7)We show by recurrence that { K,S_n } = 0 for all n ≥ 1. First, we note that as { K,H₀ } = 0, we also have { K,H₁ } = 0. As K is independent of λ_k, we also have for all G$\mathrm{⟨}\mathrm{\left\{}\mathit{K,G}\mathrm{\right\}}{\mathrm{⟩}}_{\mathit{\lambda }}\mathrm{=}\mathrm{\{}\mathit{K,}\mathrm{⟨}\mathit{G}{\mathrm{⟩}}_{\mathit{\lambda }}\mathrm{\}}\mathit{.}$(B.8)This can be seen by formal expansion of G in a Fourier series G = ∑ g_ke^{i ⟨ k,λ ⟩}. We have thus ⟨ G ⟩ _λ = g₀. Let us assume now that { K,S_k } = 0 for all k ≤ n. As ℛ_{n + 1} ∈ ℒ(H₀,H₁,S₁,...,S_n), we also have { K,ℛ_{n + 1} } = 0, and from Eq. (B.8) { K, ⟨ ℛ_{n + 1} ⟩ _λ } = 0 and thus $\mathrm{\left\{}\mathit{K,}\mathrm{\right\{}{\mathit{H}}_{\mathrm{0}}\mathit{,}{\mathit{S}}_{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{\left\}}\mathrm{\right\}}\mathrm{=}\mathrm{0.}$(B.9)Using the Jacobi identity, as { F,H₀ } = 0, we have $\mathrm{\left\{}{\mathit{H}}_{\mathrm{0}}\mathit{,}\mathrm{\right\{}\mathit{K,}{\mathit{S}}_{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{\left\}}\mathrm{\right\}}\mathrm{=}\mathrm{0.}$(B.10)The solution of the homological Eq. (B.7) is unique up to a term independent of λ. But as ⟨ { K,S_{n + 1} } ⟩ _λ = 0, then the only possible solution for Eq. (B.10) is $\mathrm{\left\{}\mathit{K,}{\mathit{S}}_{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{\right\}}\mathrm{=}\mathrm{0.}$(B.11)In the same way, as ${\mathit{H}}_{\mathrm{1}}^{\mathrm{\prime }}\mathrm{=}\mathrm{⟨}{\mathit{H}}_{\mathrm{1}}{\mathrm{⟩}}_{\mathit{\lambda }}$, we have $\mathrm{\left\{}\mathit{K,}{\mathit{H}}_{\mathrm{1}}^{\mathrm{\prime }}\mathrm{\right\}}\mathrm{=}\mathrm{⟨}\mathrm{\left\{}\mathit{K,}{\mathit{H}}_{\mathrm{1}}\mathrm{\right\}}{\mathrm{⟩}}_{\mathit{\lambda }}\mathrm{=}\mathrm{0}$. Thus { K, { H₀,S₁ } } = 0, by the Jacobi identity, { H₀, { K,S₁ } } = 0, and as previously, { K,S₁ } = 0. Our recurrence is thus complete and it follows immediately that { K,H′ } = 0.

Appendix C: Special values of C${\mathrm{c}}_{}\mathrm{(}\mathit{\alpha ,}\gtrsim \mathrm{)}$

This appendix provides the detailed computations and proofs of the results of Sect. 3.2.

Appendix C.1: Asymptotic value of C${{c}}_{}\mathrm{(}\mathit{\alpha ,}\gtrsim \mathrm{)}$ for ≳ → 0

We have shown that e_c(α, ≳) is a differentiable function of ≳, which is monotonic (41) and bounded (e_c(a, ≳ ) ∈ [ 0,1 ]). The limit ${\mathit{e}}_{\mathit{c}}\mathrm{\left(}\mathit{\alpha ,}\mathrm{0}\mathrm{\right)}\mathrm{=}{}^{\underset{\gtrsim \mathrm{\to }\mathrm{0}}{\lim }}\mathit{e}_{\mathit{c}}\mathrm{(}\mathit{\alpha ,}\gtrsim \mathrm{)}$ exists for all α ∈ ] 0,1 ]. If e_c(α, ≳) is a solution of equation(35), it will also be a solution of the following cubic equation (in e), which is directly obtained from (35) by squaring, multiplication, and simplification by α(1 + e):$\begin{array}{ccc}\mathit{K}\mathrm{\left(}\mathit{e,\alpha ,}\gtrsim \mathrm{\right)}& \mathrm{=}& \mathit{\alpha }\mathrm{(}{\gtrsim }^{\mathrm{2}}\mathrm{-}\mathit{\alpha }\mathrm{)}{\mathit{e}}^{\mathrm{3}}\mathrm{-}\mathrm{(}\mathrm{2}\mathrm{-}\mathit{\alpha }\mathrm{)}\mathrm{(}{\gtrsim }^{\mathrm{2}}\mathrm{-}\mathit{\alpha }\mathrm{)}{\mathit{e}}^{\mathrm{2}}\\ & & \end{array}$(C.1)As e_c(α, ≳), is a continuous function of ≳, when ≳ → 0 the limits e_0c(α) will satisfy the limit equation ${\mathit{K}}_{\mathrm{0}}\mathrm{\left(}\mathit{e,\alpha }\mathrm{\right)}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\alpha }\mathrm{-}\mathit{\alpha e}{\mathrm{)}}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{e}\mathrm{)}\mathrm{=}\mathrm{0}$(C.2)with solutions e₀ = 1 /α − 1 and e₁ = 1. Depending on α, several cases are treated:

α< 1 / 2:

we have then e₀> 1, and the only possibility for e_c(α,0) is e₁ = 1; as it is the only root of Eq. (C.2) which belongs to [ 0,1 ] , we have thus $\underset{\gtrsim \mathrm{\to }\mathrm{0}}{\lim }{\mathit{e}}_{\mathit{c}}\mathrm{(}\mathit{\alpha ,}\gtrsim \mathrm{)}\mathrm{=}\mathrm{1.}$(C.3)We have then $\begin{array}{ccc}\underset{\gtrsim \mathrm{\to }\mathrm{0}}{\lim }{\mathit{e}}_{\mathit{c}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\alpha ,}\gtrsim \mathrm{\right)}& \mathrm{=}& \mathrm{1}\mathrm{-}\mathrm{2}\mathit{\alpha }\mathrm{;}\\ \underset{\gtrsim \mathrm{\to }\mathrm{0}}{\lim }{\mathit{C}}_{\mathit{c}}\mathrm{\left(}\mathit{\alpha ,}\gtrsim \mathrm{\right)}& \mathrm{=}& \mathrm{1}\mathrm{-}\mathrm{2}\sqrt{\mathit{\alpha }\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\alpha }\mathrm{)}}\mathit{.}\end{array}$(C.4)In order to study the behaviour of e_c(α, ≳) in the vicinity of ≳ = 0, we can differentiate K(e_c(α, ≳ ), ≳ ) = 0 twice, which gives $\frac{\mathrm{d}{\mathit{e}}_{\mathit{c}}}{\mathrm{d}\gtrsim }\mathrm{\left(}\mathit{\alpha ,}\mathrm{0}\mathrm{\right)}\mathrm{=}\mathrm{0}\mathrm{;} \frac{{\mathit{d}}^{\mathrm{2}}{\mathit{e}}_{\mathit{c}}}{\mathit{d}{\gtrsim }^{\mathrm{2}}}\mathrm{\left(}\mathit{\alpha ,}\mathrm{0}\mathrm{\right)}\mathrm{=}\frac{\mathrm{4}\mathrm{(}\mathit{\alpha }\mathrm{-}\mathrm{1}\mathrm{)}}{\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{2}\mathit{\alpha }{\mathrm{)}}^{\mathrm{2}}}\mathit{<}\mathrm{0}\mathit{,}$(C.5)thus e_c(α, ≳) decreases with respect to ≳ at ≳ = 0.

The second order development of C_c gives ${\mathit{C}}_{\mathit{c}}\mathrm{(}\mathit{\alpha ,}\gtrsim \mathrm{)}\mathrm{=}\mathrm{1}\mathrm{-}\mathrm{2}\sqrt{\mathit{\alpha }\mathrm{(}\mathrm{1}\mathrm{-}\mathit{a}\mathrm{)}}\mathrm{+}\gtrsim \sqrt{\mathit{\alpha }}\mathrm{-}{\gtrsim }^{\mathrm{2}}\sqrt{\frac{\mathrm{1}}{\mathit{\alpha }}\mathrm{-}\mathrm{1}}\mathrm{+}\mathrm{O}\mathrm{\left(}{\gtrsim }^{\mathrm{3}}\mathrm{\right)}\mathit{.}$(C.6)

α> 1 / 2:

in this case, e₀< 1. As e_c(α, ≳ ) ∈ ] 0,e₀ [, we have e_c(α,0) ∈ [ 0,e₀ ], which gives as the unique possibility $\underset{\gtrsim \mathrm{\to }\mathrm{0}}{\lim }{\mathit{e}}_{\mathit{c}}\mathrm{(}\mathit{\alpha ,}\gtrsim \mathrm{)}\mathrm{=}\mathit{e}0$(C.7)with $\underset{\gtrsim \mathrm{\to }\mathrm{0}}{\lim }{\mathit{e}}_{\mathit{c}}^{\mathrm{\prime }}\mathrm{(}\mathit{\alpha ,}\gtrsim \mathrm{)}\mathrm{=}\mathrm{0}\mathrm{;} \underset{\gtrsim \mathrm{\to }\mathrm{0}}{\lim }{\mathit{C}}_{\mathit{c}}\mathrm{(}\mathit{\alpha ,}\gtrsim \mathrm{)}\mathrm{=}\mathrm{0.}$(C.8)By setting ≳ = 0 in Eq. (41), we also obtain $\frac{\mathrm{d}{\mathit{e}}_{\mathit{c}}}{\mathrm{d}\gtrsim }\mathrm{\left(}\mathit{\alpha ,}\mathrm{0}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{{\mathit{e}}_{\mathrm{0}}}{{\mathit{\alpha }}^{\mathrm{3}\mathit{/}\mathrm{2}}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{0}}^{\mathrm{2}}}}\mathit{<}\mathrm{0.}$(C.9)The development of C_c gives ${\mathit{C}}_{\mathit{c}}\mathrm{(}\mathit{\alpha ,}\gtrsim \mathrm{)}\mathrm{=}\gtrsim \left(\sqrt{\mathit{\alpha }}\mathrm{-}\sqrt{\mathrm{2}\mathrm{-}\frac{\mathrm{1}}{\mathit{\alpha }}}\right)\mathrm{-}\frac{{\gtrsim }^{\mathrm{2}}}{\mathrm{2}\mathit{\alpha }\mathrm{-}\mathrm{1}}\frac{\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\alpha }{\mathrm{)}}^{\mathrm{2}}}{\mathit{\alpha }}\mathrm{+}\mathrm{O}\mathrm{\left(}{\gtrsim }^{\mathrm{3}}\mathrm{\right)}\mathit{.}$(C.10)

α = 1 / 2:

in this case, e₀ = 1, and the only solution is $\underset{\gtrsim \mathrm{\to }\mathrm{0}}{\lim }{\mathit{e}}_{\mathit{c}}\mathrm{(}\mathit{\alpha ,}\gtrsim \mathrm{)}\mathrm{=}\mathrm{1}$(C.11)and $\underset{\gtrsim \mathrm{\to }\mathrm{0}}{\lim }{\mathit{e}}_{\mathit{c}}^{\mathrm{\prime }}\mathrm{(}\mathit{\alpha ,}\gtrsim \mathrm{)}\mathrm{=}\mathrm{0}\mathrm{;} \underset{\gtrsim \mathrm{\to }\mathrm{0}}{\lim }{\mathit{C}}_{\mathit{c}}\mathrm{(}\mathit{\alpha ,}\gtrsim \mathrm{)}\mathrm{=}\mathrm{0.}$(C.12)Moreover, equation (C.1) becomes ≳ ²2e²(3 − e) = (1 − e)³. We obtain thus the asymptotic value for e_c(1 / 2, ≳) when ≳ → 0 as $\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathit{c}}\left(\frac{\mathrm{1}}{\mathrm{2}}\mathit{,}\gtrsim \right)\hspace{0.17em}\hspace{0.17em}\mathrm{~}\hspace{0.17em}\hspace{0.17em}\mathrm{(}\mathrm{4}\gtrsim {\mathrm{)}}^{\mathrm{2}\mathit{/}\mathrm{3}}\mathit{.}$(C.13)For α = 1 / 2, the development of C_c in γ contains non-polynomial terms in ≳ giving ${\mathit{C}}_{\mathit{c}}\mathrm{(}\mathrm{1}\mathit{/}\mathrm{2}\mathit{,}\gtrsim \mathrm{)}\mathrm{=}\frac{\gtrsim }{\sqrt{\mathrm{2}}}\mathrm{-}{\mathrm{2}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{3}}{\gtrsim }^{\mathrm{4}\mathit{/}\mathrm{3}}\mathrm{-}{\mathrm{2}}^{\mathrm{-}\mathrm{4}\mathit{/}\mathrm{3}}{\gtrsim }^{\mathrm{2}}\mathrm{+}\mathrm{O}\mathrm{\left(}{\gtrsim }^{\mathrm{8}\mathit{/}\mathrm{3}}\mathrm{\right)}\mathit{.}$(C.14)

Appendix C.2: Asymptotic value of C${{c}}_{}\mathrm{(}\mathit{\alpha ,}\gtrsim \mathrm{)}$ for ≳ → + ∞

This case is more simple. If e_c(α, ≳) is a solution of Eq. (35), then it will also be a solution of Eq. (C.1), and thus of $\frac{\mathit{K}\mathrm{(}\mathit{e,\alpha ,}\gtrsim \mathrm{)}}{{\gtrsim }^{\mathrm{2}}}\mathrm{=}\mathrm{0.}$(C.15)As e_c(α, ≳) is monotonic and bounded, it has a limit when ≳ → + ∞, which will verify the limit equation (C.15), when ≳ → + ∞, that is ${\mathit{K}}_{\mathrm{\infty }}\mathrm{\left(}\mathit{\alpha ,e}\mathrm{\right)}\mathrm{=}{\mathit{e}}^{\mathrm{2}}\mathrm{(}\mathrm{2}\mathrm{-}\mathit{\alpha }\mathrm{-}\mathit{\alpha e}\mathrm{)}\mathrm{=}\mathrm{0.}$(C.16)As 0 <α< 1, the only solution is e = 0, and thus $\underset{\gtrsim \mathrm{\to }\mathrm{+}\mathrm{\infty }}{\lim }{\mathit{e}}_{\mathit{c}}\mathrm{(}\mathit{\alpha ,}\gtrsim \mathrm{)}\mathrm{=}\mathrm{0}$(C.17)and $\underset{\gtrsim \mathrm{\to }\mathrm{+}\mathrm{\infty }}{\lim }{\mathit{e}}_{\mathit{c}}^{\mathrm{\prime }}\mathrm{(}\mathit{\alpha ,}\gtrsim \mathrm{)}\mathrm{=}\mathrm{1}\mathrm{-}\mathit{\alpha }\mathrm{;} \underset{\gtrsim \mathrm{\to }\mathrm{+}\mathrm{\infty }}{\lim }{\mathit{C}}_{\mathit{c}}\mathrm{(}\mathit{\alpha ,}\gtrsim \mathrm{)}\mathrm{=}\mathrm{1}\mathrm{-}\sqrt{\mathit{\alpha }\mathrm{(}\mathrm{2}\mathrm{-}\mathit{\alpha }\mathrm{)}}$(C.18)and more precisely $\begin{array}{ccc}\mathit{C}\mathrm{\left(}\mathit{\alpha ,}\gtrsim \mathrm{\right)}& \mathrm{=}& \mathrm{1}\mathrm{-}\sqrt{\mathit{\alpha }\mathrm{(}\mathrm{2}\mathrm{-}\mathit{\alpha }\mathrm{)}}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}\gtrsim }\frac{\sqrt{\mathit{\alpha }}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\alpha }{\mathrm{)}}^{\mathrm{2}}}{\mathrm{2}\mathrm{-}\mathit{\alpha }}\\ & & \mathrm{+}\frac{\mathrm{1}}{{\gtrsim }^{\mathrm{2}}}\frac{\sqrt{\mathit{\alpha }\mathrm{(}\mathrm{2}\mathrm{-}\mathit{\alpha }\mathrm{)}}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\alpha }{\mathrm{)}}^{\mathrm{2}}}{\mathrm{2}\mathrm{(}\mathrm{2}\mathrm{-}\mathit{\alpha }{\mathrm{)}}^{\mathrm{3}}}\mathrm{+}\mathrm{O}\left(\frac{\mathrm{1}}{{\gtrsim }^{\mathrm{3}}}\right)\mathrm{·}\end{array}$(C.19)

Appendix C.3: Study of C${{c}}_{}\mathrm{(}\mathit{\alpha ,}\gtrsim \mathrm{)}$ for ≳ = 1 and $\gtrsim \mathrm{=}\sqrt{\mathit{\alpha }}$

For ≳ = 1 or $\gtrsim \mathrm{=}\sqrt{\mathit{\alpha }}$, we can also obtain simple expressions for C_c(α, ≳). Indeed, If ≳ = 1, K(e,α, ≳) factorises in (1 − α)(1 + e)(αe² − 2e + 1 − α) and we have a single solution for e_c in the interval [ 0,1 ], $\begin{array}{ccc}{\mathit{e}}_{\mathit{c}}\mathrm{\left(}\mathit{\alpha ,}\mathrm{1}\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{1}\mathrm{-}\sqrt{\mathrm{1}\mathrm{-}\mathit{\alpha }\mathrm{+}{\mathit{\alpha }}^{\mathrm{2}}}}{\mathit{\alpha }}\mathrm{;}\\ {\mathit{e}}_{\mathit{c}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\alpha ,}\mathrm{1}\mathrm{\right)}& \mathrm{=}& \sqrt{\mathrm{1}\mathrm{-}\mathit{\alpha }\mathrm{+}{\mathit{\alpha }}^{\mathrm{2}}}\mathrm{-}\mathit{\alpha }\mathrm{;}\end{array}$(C.20)and $\begin{array}{ccc}{\mathit{C}}_{\mathit{c}}\mathrm{\left(}\mathit{\alpha ,}\mathrm{1}\mathrm{\right)}& \mathrm{=}& \sqrt{\mathit{\alpha }}\mathrm{-}\frac{\sqrt{\mathit{\alpha }\mathrm{-}\mathrm{2}\mathrm{+}\mathrm{2}\sqrt{\mathrm{1}\mathrm{-}\mathit{\alpha }\mathrm{+}{\mathit{\alpha }}^{\mathrm{2}}}}}{\sqrt{\mathit{\alpha }}}\\ & & \end{array}$(C.21)For $\gtrsim \mathrm{=}\sqrt{\mathit{\alpha }}$, the cubic Eq. (C.1) reduces to $\mathit{K}\mathrm{\left(}\mathit{\alpha ,}\sqrt{\mathit{\alpha }}\mathrm{\right)}\mathrm{=}\mathrm{-}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\alpha }}^{\mathrm{2}}\mathrm{)}\mathit{e}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\alpha }{\mathrm{)}}^{\mathrm{2}}\mathrm{=}\mathrm{0}$(C.22)with the single solution ${\mathit{e}}_{\mathit{c}}\mathrm{\left(}\mathit{\alpha ,}\sqrt{\mathit{\alpha }}\mathrm{\right)}\mathrm{=}{\mathit{e}}_{\mathit{c}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\alpha ,}\sqrt{\mathit{\alpha }}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}\mathrm{-}\mathit{\alpha }}{\mathrm{1}\mathrm{+}\mathit{\alpha }}$(C.23)and ${\mathit{C}}_{\mathit{c}}\mathrm{\left(}\mathit{\alpha ,}\sqrt{\mathit{\alpha }}\mathrm{\right)}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{-}\sqrt{\mathit{\alpha }}{\mathrm{)}}^{\mathrm{2}}\mathit{.}$(C.24)