# AMD-stability and the classification of planetary systems

J. Laskar, A. Petit

► To cite this version:

J. Laskar, A. Petit. AMD-stability and the classification of planetary systems. Astronomy & Astrophysics - A&A, 2017, 605, pp.A72. 10.1051/0004-6361/201630022. hal-02178940

HAL Id: hal-02178940

https://hal.science/hal-02178940v1

Submitted on 21 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# AMD-stability and the classification of planetary systems*

J. Laskar and A. C. Petit

ASD/IMCCE, CNRS-UMR 8028, Observatoire de Paris, PSL, UPMC, 77 Avenue Denfert-Rochereau, 75014 Paris, France
e-mail: laskar@imcce.fr

Received 7 November 2016 / Accepted 23 January 2017

## ABSTRACT

We present here in full detail the evolution of the angular momentum deficit (AMD) during collisions as it was described in Laskar (2000, Phys. Rev. Lett., 84, 3240). Since then, the AMD has been revealed to be a key parameter for the understanding of the outcome of planetary formation models. We define here the AMD-stability criterion that can be easily verified on a newly discovered planetary system. We show how AMD-stability can be used to establish a classification of the multiplanet systems in order to exhibit the planetary systems that are long-term stable because they are AMD-stable, and those that are AMD-unstable which then require some additional dynamical studies to conclude on their stability. The AMD-stability classification is applied to the 131 multiplanet systems from The Extrasolar Planet Encyclopaedia database for which the orbital elements are sufficiently well known.

**Key words.** chaos – celestial mechanics – planets and satellites: dynamical evolution and stability – planets and satellites: formation – planets and satellites: general

## 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 $\Delta_{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

<footnote>* The AMD-stability coefficients of selected planetary systems are available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc.u-strasbg.fr/viz-bin/qcat?J/A+A/605/A72</footnote>
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_0, P_1, \dots, P_n$ be $n+1$ bodies of masses $m_0, m_1, \dots, m_n$ in gravitational interaction, and let $O$ be their centre of mass. For every body $P_i$, we denote $\mathbf{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 $(\mathbf{u}_i, \tilde{\mathbf{u}}_i = m_i \dot{\mathbf{u}}_i)_{i=0,n}$ with Hamiltonian

$$H = \frac{1}{2} \sum_{i=0}^{n} \frac{\|\tilde{\mathbf{u}}_i\|^2}{m_i} - G \sum_{0 \le i < j} \frac{m_i m_j}{\Delta_{ij}}, \quad (1)$$

where $\Delta_{ij} = \|\mathbf{u}_i - \mathbf{u}_j\|$, and $G$ is the constant of gravitation. The reduction of the centre of mass is achieved by using the canonical heliocentric variables of Poincaré ($\mathbf{r}_i, \tilde{\mathbf{r}}_i$) (Poincaré 1905; Laskar & Robutel 1995), defined as

$$\mathbf{r}_0 = \mathbf{u}_0; \quad \mathbf{r}_i = \mathbf{u}_i - \mathbf{u}_0 \quad \text{for } i \neq 0 \quad (2)$$

$$\tilde{\mathbf{r}}_0 = \tilde{\mathbf{u}}_0 + \tilde{\mathbf{u}}_1 + \dots + \tilde{\mathbf{u}}_n; \quad \tilde{\mathbf{r}}_i = \tilde{\mathbf{u}}_i \quad \text{for } i \neq 0. \quad (3)$$

This Hamiltonian can then be split into an integrable part, $H_0$, and a perturbation, $H_1$,

$$H = H_0 + H_1 \quad (4)$$

with

$$H_0 = \frac{1}{2} \sum_{i=1}^{n} \frac{\|\tilde{\mathbf{r}}_i\|^2}{m_i} - G \sum_{i=1}^{n} \frac{m_0 m_i}{r_i} \quad (5)$$

and

$$H_1 = \frac{1}{2} \frac{\|\tilde{\mathbf{u}}_0\|^2}{m_0} - G \sum_{1 \le i < j} \frac{m_i m_j}{\Delta_{ij}} \quad (6)$$

The integrable part, $H_0$, 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_0$. A set of adapted variables for $H_0$ will thus be given by the elliptical elements, ($a_k, e_k, i_k, \lambda_k, \omega_k, \Omega_k$), where $a_k$ is the semi-major axis, $e_k$ the eccentricity, $i_k$ the inclination, $\lambda_k$ the mean longitude, $\omega_k$ the longitude of the perihelion, and $\Omega_k$ the longitude of the node. They are defined as the elliptical elements associated to the Hamiltonian

$$H_{0k} = \frac{1}{2} \frac{\|\tilde{\mathbf{r}}_k\|^2}{m_k} - \mu \frac{m_k}{r_k} \quad (7)$$

with $\mu = Gm_0$.

### 2.2. Angular momentum deficit (AMD)

The total linear momentum of the system is null in the barycentric reference frame

$$L = \sum_{i=0}^{n} m_i \dot{\mathbf{u}}_i = \sum_{i=0}^{n} \tilde{\mathbf{u}}_i = \tilde{\mathbf{r}}_0 = 0. \quad (8)$$

Let $G$ be the total angular momentum. Its expression is not changed by the linear symplectic change of variable $(\mathbf{u}, \tilde{\mathbf{u}}) \rightarrow (\mathbf{r}, \tilde{\mathbf{r}})$. We have thus

$$G = \sum_{i=0}^{n} \mathbf{u}_i \wedge \tilde{\mathbf{u}}_i = \sum_{i=1}^{n} \mathbf{r}_i \wedge \tilde{\mathbf{r}}_i. \quad (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_0$. In particular, if the angular momentum direction is assumed to be the axis $z$, the norm of the angular momentum is

$$G = \sum_{k=1}^{n} \Lambda_k \sqrt{1 - e_k^2} \cos i_k, \quad (10)$$

where $\Lambda_k = m_k \sqrt{\mu a_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)

$$C = \sum_{k=1}^{n} \Lambda_k \left( 1 - \sqrt{1 - e_k^2} \cos i_k \right). \quad (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_1$ and $m_2$, and orbits $O_1, O_2$ occurs, forming a new body ($m_3, O_3$). During this collision we assume that the other bodies are not affected. The mass is conserved

$$m_3 = m_1 + m_2, \quad (12)$$

and the linear momentum in the barycentric reference frame is conserved so $\tilde{\mathbf{u}}_3 = \tilde{\mathbf{u}}_1 + \tilde{\mathbf{u}}_2$, and also

$$\tilde{\mathbf{r}}_3 = \tilde{\mathbf{r}}_1 + \tilde{\mathbf{r}}_2. \quad (13)$$

On the other hand, $\mathbf{r}_3 = \mathbf{r}_1 = \mathbf{r}_2$ at the time of the collision, so the angular momentum is also conserved

$$\mathbf{r}_3 \wedge \tilde{\mathbf{r}}_3 = \mathbf{r}_1 \wedge \tilde{\mathbf{r}}_1 + \mathbf{r}_2 \wedge \tilde{\mathbf{r}}_2. \quad (14)$$

It should be noted that the transformation of the orbits ($m_1, O_1$) + ($m_2, O_2$) $\rightarrow$ ($m_3, O_3$) 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 ($m_1, O_1$) and ($m_2, O_2$) 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

$$h_i = \frac{1}{2} \frac{\|\tilde{\mathbf{r}}_i\|^2}{m_i} - \mu \frac{m_i}{r_i} = -\mu \frac{m_i}{2a_i}. \quad (15)$$
At collision, we have $\mathbf{r}_1 = \mathbf{r}_2 = \mathbf{r}_3 = \mathbf{r}$; we have thus the conservation of the potential energy

$$-\mu \frac{m_3}{r_3} = -\mu \frac{m_1 + m_2}{r} = -\mu \frac{m_1}{r_1} - \mu \frac{m_2}{r_2}. \quad (16)$$

The change of Keplerian energy is thus given by the change of kinetic energy

$$\delta h = h_3 - h_1 - h_2 = \frac{1 \|\tilde{\mathbf{r}}_3\|^2}{2 m_3} - \frac{1 \|\tilde{\mathbf{r}}_1\|^2}{2 m_1} - \frac{1 \|\tilde{\mathbf{r}}_2\|^2}{2 m_2}; \quad (17)$$

that is, with Eqs. (12), (13),

$$\begin{aligned} 2m_1 m_2 m_3 \delta h &= m_1 m_2 (\tilde{\mathbf{r}}_1 + \tilde{\mathbf{r}}_2)^2 - m_2 (m_1 + m_2) \tilde{\mathbf{r}}_1^2 - m_1 (m_1 + m_2) \tilde{\mathbf{r}}_2^2 \\ &= -m_2^2 \tilde{\mathbf{r}}_1^2 - m_1^2 \tilde{\mathbf{r}}_2^2 + 2m_1 m_2 \tilde{\mathbf{r}}_1 \cdot \tilde{\mathbf{r}}_2 \\ &= -(m_2 \tilde{\mathbf{r}}_1 - m_1 \tilde{\mathbf{r}}_2)^2 \\ &\le 0. \end{aligned} \quad (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 $m_2 \tilde{\mathbf{r}}_1 = m_1 \tilde{\mathbf{r}}_2$, that is, as $m_1 m_2 \neq 0$, when $\dot{\mathbf{u}}_1 = \dot{\mathbf{u}}_2$. As an immediate consequence of the decrease of energy during the collision, we have

$$\frac{1}{a_3} \ge \frac{m}{a_1} + \frac{1-m}{a_2} \quad (19)$$

with

$$m = \frac{m_1}{m_1 + m_2} = \frac{m_1}{m_3}. \quad (20)$$

### 2.3.2. AMD evolution during collision

Let $f(x) = 1/\sqrt{x}$. As $f'(x) < 0$ and $f''(x) > 0$, we have, as $f$ is decreasing and convex,

$$f\left(\frac{1}{a_3}\right) \le f\left(\frac{m}{a_1} + \frac{1-m}{a_2}\right) \le m f\left(\frac{1}{a_1}\right) + (1-m) f\left(\frac{1}{a_2}\right) \quad (21)$$

thus

$$m_3 \sqrt{a_3} \le m_1 \sqrt{a_1} + m_2 \sqrt{a_2}. \quad (22)$$

During the collision, the angular momentum is conserved Eq. (14), and so is the conservation of its normal component, that is

$$\begin{aligned} m_3 \sqrt{\mu a_3} \sqrt{1 - e_3^2} \cos i_3 &= m_1 \sqrt{\mu a_1} \sqrt{1 - e_1^2} \cos i_1 + m_2 \sqrt{\mu a_2} \sqrt{1 - e_2^2} \cos i_2. \end{aligned} \quad (23)$$

We deduce that in all circumstances we have a decrease of the angular momentum deficit during the collision, that is

$$C_3 \le C_1 + C_2 \quad (24)$$

with

$$C_k = m_k \sqrt{\mu a_k} (1 - \sqrt{1 - e_k^2} \cos i_k) \quad (k = 1, 3). \quad (25)$$

The equality can hold in Eq. (24) only if $m_1 = 0$, $m_2 = 0$, or $a_1 = a_2$ and $\dot{\mathbf{u}}_1 = \dot{\mathbf{u}}_2$, 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 $\mathbf{r}_1 = \mathbf{r}_2$).

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<fnref n="1" /> 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

$$a(1 + e) = a'(1 - e') \quad (26)$$

$$m \sqrt{\mu a} (1 - \sqrt{1 - e^2}) + m' \sqrt{\mu a'} (1 - \sqrt{1 - e'^2}) = C, \quad (27)$$

where $(m, a, e)$ are the parameters of the inner orbit and $(m', a', e')$ those of the outer orbit ($a \le a'$).

### 3.1. Collisional condition

We assume that $a, a', m, m'$ are non-zero. Denoting $\alpha = a/a'$, $\gamma = m/m'$, the system in Eq. (27) becomes

$$\mathcal{D}(e, e') = \alpha e + e' - 1 + \alpha = 0 \quad (28)$$

$$\begin{aligned} \mathcal{C}(e, e') &= \gamma \sqrt{\alpha} (1 - \sqrt{1 - e^2}) + (1 - \sqrt{1 - e'^2}) \\ &= C/\Lambda' \end{aligned} \quad (29)$$

with $\Lambda' = m' \sqrt{\mu a'}$, and where $\mathcal{C}(e, e') = C/\Lambda'$ is called the relative AMD. As $e$ and $e'$ are planetary eccentricities, we also have

$$0 \le e \le 1; \quad 0 \le e' \le 1. \quad (30)$$

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 $\mathcal{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 $e'_0 = 1 - \alpha < 1$, and the intersection of the collision line with the axis $e' = 0$ is obtained for $e_0 = 1/\alpha - 1$. This value can be greater or smaller than 1 depending on the value of $\alpha$. We have thus the different cases (Fig. 1)

$$(a) \quad \alpha < \frac{1}{2} \iff e_0 > 1 \quad \text{and} \quad e \le 1 \quad (31)$$

$$(b) \quad \alpha > \frac{1}{2} \iff e_0 < 1 \quad \text{and} \quad e < e_0 \quad (32)$$

<footnote n="1">¹ Here practical stability means stability over a very long time compared to the expected life of the central star.</footnote>
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 1. Collision conditions for $e_0 = 1/\alpha - 1$. Case (a): $\alpha < 1/2 \Leftrightarrow e_0 > 1$. Case (b): $\alpha > 1/2 \Leftrightarrow e_0 < 1$. Collisions can only occur in the shaded region.</figcaption></figure>

and the limit case $\alpha = 1/2$, for which $e_0 = 1$. In all cases, the Lagrange multipliers condition is written

$$ \lambda \nabla \mathcal{D}(e, e') = \nabla \mathcal{C}(e, e'), \quad (33) $$

which gives

$$ \frac{\sqrt{1 - e'^2}}{e'} = \frac{\sqrt{\alpha} \sqrt{1 - e^2}}{\gamma e}. \quad (34) $$

This relation allows $e'$ to be eliminated in the collision condition (28), which becomes an equation in the single variable $e$, and parameters $(\alpha, \gamma)$:

$$ F(e, \alpha, \gamma) = \alpha e + \frac{\gamma e}{\sqrt{\alpha(1 - e^2) + \gamma^2 e^2}} - 1 + \alpha = 0. \quad (35) $$

Here $F(e, a, \gamma)$ is properly defined for $(e, a, \gamma)$ in the domain $D_{e,a,\gamma}$ defined by $e \in [0, 1]$, $\alpha \in [0, 1]$, $\gamma \in [0, +\infty[$, as in this domain, $1 - e^2 + \gamma^2 e^2 / \alpha > 0$. We also have

$$ \frac{\partial F}{\partial e}(e, \alpha, \gamma) = \alpha + \frac{\alpha \gamma}{(\alpha(1 - e^2) + \gamma^2 e^2)^{3/2}}. \quad (36) $$

Thus, $\frac{\partial F}{\partial e}(e, \alpha, \gamma) > 0$ on the domain $D_{e,a,\gamma}$ and $F(e, \alpha, \gamma)$ is strictly increasing with respect to $e$ for $e \in [0, 1]$. Moreover, as $0 < \alpha < 1$,

$$ F(0, \alpha, \gamma) = -1 + \alpha < 0; \quad F(1, \alpha, \gamma) = 2\alpha > 0; \quad (37) $$

and

$$ F(e_0, \alpha, \gamma) = \frac{\gamma e_0}{\sqrt{\alpha(1 - e_0^2) + \gamma^2 e_0^2}} > 0. \quad (38) $$

The equation of collision (35) thus always has a single solution $e_c$ in the interval $[0, \min(1, e_0)]$. This ensures that this critical value of $e$ will fulfil the condition (30). The corresponding value of the relative AMD $C_c(\alpha, \gamma) = \mathcal{C}(e_c, e'_c)$ is then obtained through (29).

## 3.2. Critical AMD $C_c(\alpha, \gamma)$

We have thus demonstrated that for a given pair of ratios of semi-major axes and masses, $(\alpha, \gamma)$, there is always a unique critical value $C_c(\alpha, \gamma)$ of the relative AMD $\mathcal{C} = C/\Lambda'$ which defines the AMD-stability. The system of two planets is AMD-stable if and only if

$$ \mathcal{C} = \frac{C}{\Lambda'} < C_c(\alpha, \gamma). \quad (39) $$

The value of the critical AMD $C_c(\alpha, \gamma)$ is obtained by computing first the critical eccentricity $e_c(\alpha, \gamma)$ which is the unique solution of the collisional Eq. (35) in the interval $[0, 1]$. The critical AMD is then $C_c(\alpha, \gamma) = \mathcal{C}(e_c, e'_c)$ (Eq. (29)), where the critical value $e'_c$ 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 $(\alpha, \gamma)$.

## 3.3. Behaviour of the critical AMD

We now analyse the general properties of the critical AMD function $C_c(\alpha, \gamma)$. As $\frac{\partial F}{\partial e}(e, \alpha, \gamma) > 0$, we can apply the implicit function theorem to the domain $D_{e,a,\gamma}$, which then ensures that the solution of the collision Eq. (35), $e_c(\alpha, \gamma)$, is a continuous function of $\gamma$ (and even analytic for $\gamma \in [0, +\infty[$). Moreover, on $D_{e,a,\gamma}$,

$$ \frac{\partial F}{\partial \gamma}(e, \alpha, \gamma) = \frac{\alpha e(1 - e^2)}{(\alpha(1 - e^2) + \gamma^2 e^2)^{3/2}} \geq 0. \quad (40) $$

We also have

$$ \frac{\partial e_c}{\partial \gamma}(\alpha, \gamma) = - \frac{e_c(1 - e_c^2)}{(\alpha(1 - e_c^2) + \gamma^2 e_c^2)^{3/2} + \gamma} \le 0, \quad (41) $$

and $e_c(\alpha, \gamma)$ is a decreasing function of $\gamma$. For any given values of the semi-major axes ratio $\alpha$, and masses $\gamma$, we can thus find the critical value $C_c(\alpha, \gamma)$ which allows for a collision (39). For the critical value $C_c(\alpha, \gamma)$, a single solution corresponds to the tangency condition (Fig. 1), and this solution is obtained at the critical value $e_c(\alpha, \gamma)$ for the eccentricity of the orbit $O$. The values of the critical relative AMD $C_c(\alpha, \gamma)$ are plotted in Fig. 2
<table><caption>**Table 1.** Special values of $C_c(\alpha, \gamma)$.</caption><thead><tr><th>γ</th><th>α</th><th>e<sub>c</sub>(α, γ)</th><th>e'<sub>c</sub>(α, γ)</th><th>C<sub>c</sub>(α, γ)</th></tr></thead><tbody><tr><td>γ → 0</td><td>α < 1/2</td><td>1 - 2(1 - α)/(1 - 2α)<sup>2</sup>γ<sup>2</sup></td><td>1 - 2α + 2α(1 - α)/(1 - 2α)<sup>2</sup>γ<sup>2</sup></td><td>1 - 2√α(1 - α) + √αγ</td></tr><tr><td>γ → 0</td><td>α = 1/2</td><td>1 - (4γ)<sup>2/3</sup></td><td>2<sup>1/3</sup>γ<sup>2/3</sup></td><td>γ/√2</td></tr><tr><td>γ → 0</td><td>α > 1/2</td><td>e<sub>0</sub> - <sup>e<sub>0</sub></sup>⁄<sub>√α(2α - 1)γ</sub></td><td><sup>√αe<sub>0</sub>γ</sup>⁄<sub>√2α - 1</sub></td><td>(√α - √2 - 1/α)γ</td></tr><tr><td>γ → +∞</td><td>0 < α < 1</td><td><sup>1</sup>⁄<sub>γ</sub>(1 - α)/√2 - α</td><td>1 - α - <sup>α(1 - α)</sup>⁄<sub>γ</sub>√2 - α</td><td>1 - √α(2 - α) - <sup>√α(1 - α)<sup>2</sup></sup>⁄<sub>2 - α</sub>(1/γ)</td></tr><tr><td>1</td><td>0 ≤ α ≤ 1</td><td><sup>1 - √1 - α + α<sup>2</sup></sup>⁄<sub>α</sub></td><td>√1 - α + α<sup>2</sup> - α</td><td>1 + √α - <sup>√α - 2 + 2√1 - α + α<sup>2</sup></sup>⁄<sub>√α</sub><br/>- √α√1 - 2α + 2√1 - α + α<sup>2</sup></td></tr><tr><td>√α</td><td>0 ≤ α ≤ 1</td><td><sup>1 - α</sup>⁄<sub>1 + α</sub></td><td><sup>1 - α</sup>⁄<sub>1 + α</sub></td><td>(1 - √α)<sup>2</sup></td></tr></tbody></table>

**Notes.** The detail of the computations is provided in Appendix C.

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Values of $C_c(\alpha, \gamma)$ vs. $\alpha$ for the different $\gamma$ values for which an explicit expression of $C_c(\alpha, \gamma)$ is obtained.</figcaption></figure>

versus $\alpha$, for different values of $\gamma$. Deriving Eq. (29) with respect to $\gamma$, one obtains

$$ \frac{\partial C}{\partial \gamma} = \sqrt{\alpha} (1 - \sqrt{1 - e^2}) + \gamma \sqrt{\alpha} \frac{e}{\sqrt{1 - e^2}} \frac{\partial e}{\partial \gamma} + \frac{e'}{\sqrt{1 - e'^2}} \frac{\partial e'}{\partial \gamma}. \quad (42) $$

Using the two relations (Eqs. (28), (34)), this reduces to

$$ \frac{\partial C_c(\alpha, \gamma)}{\partial \gamma} = \sqrt{\alpha} (1 - \sqrt{1 - e_c^2}) > 0. \quad (43) $$

Thus $C_c(\alpha, \gamma)$ strictly increases with $\gamma$. In the same way

$$ \frac{\partial C_c(\alpha, \gamma)}{\partial \alpha} = \frac{\gamma}{2\sqrt{\alpha}} \left( 1 - \frac{(1 + e_c^2)}{\sqrt{1 - e_c^2}} \right) < 0 \quad (44) $$

and $C_c(\alpha, \gamma)$ decreases with $\alpha$.

Now that the general behaviour of the critical AMD $C_c(\alpha, \gamma)$ 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(\alpha, \gamma)$ for $\alpha \to 1$ can also be obtained (see Appendix C). With $\eta = 1 - \alpha$, we have

$$ C_c(\alpha, \gamma) \sim \frac{\gamma}{\gamma + 1} \frac{\eta^2}{2} + O(\eta^3). \quad (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{1}{2}m\sqrt{ae^2}$. As the mass of the particle increases, its ex-
cursion 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 plan-
etesimal with semi-major axis *a*′, such that there is no other plan-
etesimal 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*<sub>c</sub>(*α*, *γ*):

$$
C_{c1}(\alpha, \gamma) = k(\gamma) \left(\frac{\delta_a}{a}\right)^2, \qquad (46)
$$

where $\delta_a = a' - a$ and $k(\gamma) = \gamma/(2(\gamma + 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 planetes-
imals $\rho(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 planetesi-
mal 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

$$
\mathcal{C} = \frac{C}{\Lambda'} \ge C_{c1} = k \left(\frac{\delta_a}{a}\right)^2 . \qquad (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

$$
m(a') \sim m(a) \sim \rho(a) \delta_a, \quad (48)
$$

and from (47) at the limiting case we have

$$
\frac{\tilde{C}}{\delta_a \rho(a) \sqrt{a}} = k \left(\frac{\delta_a}{a}\right)^2, \qquad (49)
$$

where $\tilde{C} = C/\sqrt{\mu}$. We have thus

$$
\delta_a = \left(\frac{\tilde{C}}{k}\right)^{1/3} a^{1/2} \rho^{-1/3}, \qquad (50)
$$

and from Eq. (48)

$$
m(a) = \left( \frac{\tilde{C}}{k} \right)^{1/3} a^{1/2} \rho^{2/3}. \qquad (51)
$$

4.3. Scaling laws with initial mass distribution $\rho(a) \propto a^p$

Using the previous relations, we can compute the resulting sys-
tems for various initial mass distribution, in particular for ρ(a) =
ζa^p. From Eq. (51), we obtain for two consecutive planets

$$
\frac{m}{m'} = \alpha^{1/2} \left( \frac{\rho(a)}{\rho(a')} \right)^{2/3} = \alpha^{\frac{1}{2} + \frac{1}{3} p} \qquad (52)
$$

and from Eq. (46), as $\lim_{\alpha \to 1} \gamma = 1$, we thus have $k(\gamma) = \frac{1}{4}$ in all cases. Relation (50) can be rewritten

$$
\delta_a a^{\frac{p}{3}-\frac{1}{2}} = (4\tilde{C})^{1/3} \zeta^{-1/3} \delta n, \quad (53)
$$

where $\delta n = 1$ is the increment from planet *a* to *a*'. By integration,
this difference relation becomes for $p \neq -\frac{3}{2}$,

$$
a^{-\frac{2p+3}{6}} = a_{0}^{-\frac{2p+3}{6}} + \frac{2p+3}{6} \left( \frac{4\tilde{C}}{\zeta} \right)^{1/3} n. \quad (54)
$$

For $p = -\frac{3}{2}$, we obtain

$$
\log(a) = \log(a_0) + (4\tilde{C})^{1/3} \zeta^{-1/3} n. \quad (55)
$$

In particular, for $p=0$ (constant distribution and $\rho(a) = \zeta$), from Eq. (53), we have

$$
\sqrt{a} = \sqrt{a_0} + \left(\frac{\tilde{C}}{2\zeta}\right)^{1/3} n. \qquad (56)
$$

For the masses, from Eq. (51), we have for large *n* (or if *a*₀ is small)

$$
m(n) \sim (2\tilde{C}^2\zeta)^{1/3}n. \qquad (57)
$$

For $p = -3/2$, we find a power law similar to the Titius-Bode law for planetary distribution<fnref n="2" />. 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 intersec-
tions. When an intersection occurs, the two bodies merge into
a new one whose orbital parameters are determined by the col-
lisional Eqs. (12), (13) (see Appendix A). Between collisions,
the orbits do not evolve, apart from a diffusion of their eccen-
tricities, 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).

<footnote n="2">² The value *p* = -3/2 corresponds to a surface density proportional to *r*⁻⁵/², which is different from the -3/2 surface density exponent of the minimum mass solar nebula (Weidenschilling 1977). For the solar system, (Laskar 2000) found *p* = 0 to be the best fitting value; it corresponds to a surface density proportional to *r*⁻¹.</footnote>
<table><caption>**Table 2.** Planetary distribution corresponding to different initial mass distribution.</caption><thead><tr><th>p</th><th>a(n)</th><th>m(a)</th><th>m(n)</th></tr></thead><tbody><tr><td>p ≠ -3/2</td><td>a<sup>2p+3/6</sup> = a<sub>0</sub><sup>2p+3/6</sup> + <sup>2p+3</sup>&frasl;<sub>6</sub> (4&#x0303;C/&zeta;)<sup>1/3</sup> n</td><td>(4&#x0303;C&zeta;<sup>2</sup>)<sup>1/3</sup> a<sup>4p+3/6</sup></td><td>m(n) ~ 4&#x0303;C (4&#x0303;C/&zeta;)<sup>-1/2p+3/6</sup> (2p + 3/6) n<sup>2-3/2p+3/6</sup></td></tr><tr><td>-3/2</td><td>log(a) = log(a<sub>0</sub>) + (4&#x0303;C/&zeta;)<sup>1/3</sup> n</td><td>(4&#x0303;C&zeta;<sup>2</sup>)<sup>1/3</sup> a<sup>-1/2</sup></td><td>log(m) ~ -(&#x0303;C/2&zeta;)<sup>1/3</sup> n</td></tr><tr><td>0</td><td>&radic;a = &radic;(a<sub>0</sub> + (4&#x0303;C/&zeta;)<sup>1/3</sup> n<sup>2</sup>/2)</td><td>(4&#x0303;C&zeta;<sup>2</sup>)<sup>1/3</sup> a<sup>1/2</sup></td><td>m ~ (2&#x0303;C<sup>2</sup>&zeta;)<sup>1/3</sup> n</td></tr></tbody></table>

# 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

$$ \beta = \frac{C}{C_c} = \frac{C}{\Lambda'C_c} \qquad (58) $$

We call $\beta$ the AMD-stability coefficient. For pairs of planets, $\beta < 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<fnref n="3" />. We can also define an AMD-stability coefficient regarding the collision with the star. We define $\beta_S$, the AMD-stability coefficient of the pair formed by the star and the innermost planet. For this pair, we have $\alpha = 0$ and thus $C_c = 1$. With this simplification $\beta_S = C/\Lambda$, where $\Lambda$ 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<fnref n="4" /> (Wright et al. 2011).

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 $\alpha$ and $\gamma$, we can use equivalent quantities that are measured more precisely in observations than the masses

<figure><img src="image_5.png" /><figcaption>**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.</figcaption></figure>

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 $m \sin(i)$ for RV system. This is not a problem for the computation of $\gamma$; if we assume that the systems are close to coplanarity, then

$$ \gamma = \frac{m}{m'} = \frac{m \sin(i)}{m' \sin(i)} \approx \frac{m \sin(i)}{m' \sin(i')} \qquad (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

$$ \beta_p = C_p/(\Lambda'C_c), \qquad (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

$$ \beta = 2\beta_p. \qquad (61) $$

<footnote n="3">³ This is equivalent to require that all pairs are AMD-stable.</footnote>

<footnote n="4">⁴ http://exoplanets.org/</footnote>
<table><caption>**Table 3.** Result of the analysis split in function of the multiplicity of the system.</caption><thead><tr><th>Multiplicity</th><th>Strong stable</th><th>Weak stable</th><th>Unstable</th><th>Total</th></tr></thead><tbody><tr><td>2</td><td>42</td><td>21</td><td>34</td><td>97</td></tr><tr><td>3</td><td>4</td><td>1</td><td>17</td><td>22</td></tr><tr><td>4+</td><td>2</td><td>0</td><td>10</td><td>12</td></tr><tr><td>Total</td><td>48</td><td>22</td><td>61</td><td>131</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Probability distribution function of $\beta$ for the sample studied. The systems are grouped by multiplicity. The vertical line $\beta = 1$ marks the separation between AMD-stable and AMD-unstable pairs.</figcaption></figure>

We can also define a coplanar AMD-stability coefficient associated with the star, and similarly we set $\beta_S = 2\beta_{Sp}$. We then compute the coefficients $\beta_S$ and $\beta$ 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 $\beta$ for the considered systems.

## 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 $(\alpha, \gamma, C_c, \beta...)$ in these 10 000 cases and the associated distributions.

For masses (or $m \sin(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 $(e \cos \omega, e \sin \omega)$ are Gaussian. Since the average value of $\omega$ has no importance in the computation of the eccentricity distribution, we assume it to be 0. Therefore, $e \sin \omega$ has 0 mean. We define the distribution of $\tilde{e} = e \cos \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 $e \sin \omega$ has the same standard deviation as $e \cos \omega$, we have $e \sin \omega = \tilde{e} - \langle \tilde{e} \rangle$. The distribution of $e$ is then deduced from that of $\tilde{e}$ using

$$e = \sqrt{\tilde{e}^2 + (\tilde{e} - \langle \tilde{e} \rangle)^2}. \quad (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 $\alpha$ or $\gamma$ 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 $\alpha$,

$$\alpha = \frac{a}{a'} = \frac{a}{\langle a' \rangle} \left( 1 - \frac{\Delta a'}{\langle a' \rangle} \right), \quad (63)$$

with $\Delta a' = a' - \langle a' \rangle$.

## 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 $\beta < 1$). In addition, if the AMD-stability coefficient of the star $\beta_S < 1$ (resp. $\beta_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.
<figure><img src="image_7.png" /><figcaption>**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.</figcaption></figure>

### 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.
<figure><img src="image_8.png" /><figcaption>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 $\beta_S$).</figcaption></figure>

### 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 $H_0$, 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.

We first test $H_0$ 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 $H_0$ with a p-value of about 9%. Therefore, we cannot reject $H_0$ 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 $F$. Let us consider an interval $I_x = (x, x + \Delta x)$, if we assume $H_0$, the probability for a ratio r to be in $I_x$ is

$$ \begin{aligned} \mathcal{P}(r \in I_x | H_0) &= \int_{x}^{x+\Delta x} f(r') dr' \\ &= \mathcal{F}(x+\Delta x) - \mathcal{F}(x) = \Delta \mathcal{F}(x). \end{aligned} \quad (64) $$

Therefore, the probability that more than $k_0$ pairs out of N pairs drawn randomly from the distribution $F$ are in $I_x$ is

$$ \mathcal{P}(x, k_0 | H_0) = \sum_{k=k_0}^{N} \binom{N}{k} (1 - \Delta \mathcal{F})^{N-k} \Delta \mathcal{F}^k. \quad (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 $\Delta x$, we can compute the probability $P(x, R_u) = \mathcal{P}(x, k_u(x))$, 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 $F$. We choose $\Delta 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 \times 10^{-3}$ and $2.9 \times 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 $F$. We create 10,000 datasets $R_\alpha$ by drawing N points from $F$ and compute $P(x, R_\alpha)$ for these datasets. Then, we record the two minimum values (at least distant by more than $\Delta 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.

## 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
<figure><img src="image_9.png" /><figcaption>**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$).</figcaption></figure>

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 obtained 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.
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Circular momenta and AMD of the solar system planets. The horizontal lines correspond to the respective mean $\Lambda_k$ and $C_k$ for the inner and outer system.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Examples of systems where the outer part is AMD-stable and the inner part becomes AMD-stable if considered alone.</figcaption></figure>

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 stabilize 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

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Cumulative period ratios of the AMD-unstable systems and of all systems in the catalogue. The vertical line represents low order MMR.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Red curve: probability of observing more than $k_0$ pairs in the interval $I_x$ given a probability density $f$. Green curve: $k_0$ is the number of pairs in $I_x$ for the unstable systems $R_u$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Probability density function of the two first minimum values of $P(x, R_\alpha)$ from 10 000 samples $R_\alpha$ drawn from the $\mathcal{F}$ distribution.</figcaption></figure>

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 consider-
ation 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 consider-
ing 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 al-
lowing 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

Albouy, A. 2002, Lectures on the two-body problem (Princeton University Press)

Batygin, K., Morbidelli, A., & Holman, M. J. 2015, ApJ, 799, 120

Chambers, J. 2001, Icarus, 152, 205

Chambers, J. 2013, Icarus, 224, 43

Chambers, J., Wetherill, G., & Boss, A. 1996, Icarus, 119, 261

Deck, K. M., Payne, M., & Holman, M. J. 2013, ApJ, 774, 129

Fabrycky, D. C., Lissauer, J. J., Ragozzine, D., et al. 2014, ApJ, 790, 146

Fang, J., & Margot, J.-L. 2012, ApJ, 761, 92

Fang, J., & Margot, J.-L. 2013, ApJ, 767, 115

Figueira, P., Marmier, M., Boué, G., et al. 2012, A&A, 541, A139

Gladman, B. 1993, Icarus, 106, 247

Graner, F., & Dubrulle, B. 1994, A&A, 282, 262

Hernández-Mena, C., & Benet, L. 2011, MNRAS, 412, 95

Kokubo, E., & Genda, H. 2010, ApJ, 714, L21

Laskar, J. 1994, A&A, 287, L9

Laskar, J. 1997, A&A, 317, L75

Laskar, J. 2000, Phys. Rev. Lett., 84, 3240

Laskar, J. 2008, Icarus, 196, 1

Laskar, J., & Gastineau, M. 2009, Nature, 459, 817

Laskar, J., & Robutel, P. 1995, Celest. Mech. Dyn. Astron., 62, 193

Lehmann, E. L., & Romano, J. P. 2006, Testing statistical hypotheses (Springer Science & Business Media)

Lissauer, J. J., Ragozzine, D., Fabrycky, D. C., et al. 2011, ApJS, 197, 8

Marchal, C., & Bozis, G. 1982, Celest. Mech., 26, 311

Mayor, M., Marmier, M., Lovis, C., et al. 2011, A&A, submitted
[arXiv:1109.2497]

Moorhead, A. V., Ford, E. B., Morehead, R. C., et al. 2011, ApJS, 197, 1

Morbidelli, A., Lunine, J., O'Brien, D., Raymond, S., & Walsh, K. 2012, Annu. Rev. Earth Planet. Sci., 40, 251

Nieto, M. M. 1972, The Titius-Bode law of planetary distances: its history and theory (Pergamon Press), 161

Poincaré, H. 1905, Leçons de mécanique céleste, Tome I (Paris: Gauthier-Villars)

Pu, B., & Wu, Y. 2015, ApJ, 807, 44

Ramos, X. S., Correa-Otto, J. A., & Beaugé, C. 2015, Celest. Mech. Dyn. Astron., 123, 453

Raymond, S. N., Quinn, T., & Lunine, J. I. 2006, Icarus, 183, 265

Schneider, J., Dedieu, C., Le Sidaner, P., Savalle, R., & Zolotukhin, I. 2011, A&A, 532, A79

Shabram, M., Demory, B.-O., Cisewski, J., Ford, E. B., & Rogers, L. 2016, ApJ, 820, 93

Smith, A. W., & Lissauer, J. J. 2009, Icarus, 201, 381

Tremaine, S. 2015, ApJ, 807, 157

Weidenschilling, S. J. 1977, Astrophys. Space Sci., 51, 153

Winn, J. N., & Fabrycky, D. C. 2015, Ann. Rev. Astron. Astrophys., 53, 409

Wisdom, J. 1980, AJ, 85, 1122

Wright, J. T., Fakhouri, O., Marcy, G. W., et al. 2011, PASP, 123, 412

Xie, J.-W., Dong, S., Zhu, Z., et al. 2016, Proceedings of the National Academy of Sciences, 113, 11431
## Appendix A: Intersection of planar orbits

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 $(\mu, r, \dot{r})$ and let $G = r \wedge \dot{r}$ be the angular momentum per unit of mass. The Laplace vector

$$ \mathbf{P} = \frac{\dot{r} \wedge G}{\mu} - \frac{r}{r} \qquad (A.1) $$

is an integral of the motion with coordinates $(e \cos \omega, e \sin \omega)$. One has

$$ \mathbf{P} \cdot \mathbf{r} = \frac{G^2}{\mu} - r = p - r, \qquad (A.2) $$

where $p = a(1 - e^2)$ 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

$$ r^2 = x^2 + y^2 \qquad (A.3) $$

with the plane defined by Eq. (A.2), that is

$$ x(e \cos \omega) + y(e \sin \omega) + r = p. \qquad (A.4) $$

If we consider now two orbits $O_1$ and $O_2$. Their intersection is easily obtained as the intersection of the line of intersection of the two planes

$$ \begin{aligned} & x(e_1 \cos \omega_1) + y(e_1 \sin \omega_1) + r = p_1 \\ & x(e_2 \cos \omega_2) + y(e_2 \sin \omega_2) + r = p_2 \end{aligned} \qquad (A.5) $$

with the cone of equation $r^2 = x^2 + y^2$. Depending on the initial conditions, if $O_1$ and $O_2$ are distinct, we will get either 0, 1, or 2 solutions.

<figure><img src="image_15.png" /><figcaption>Fig. A.1. Elliptical orbit, as the intersection of the cone $r^2 = x^2 + y^2$ and the plane $\mathbf{P} \cdot \mathbf{r} + r = p$.</figcaption></figure>

## 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 $\lambda_k$. Let

$$ H = H_0(\Lambda) + \varepsilon H_1(\Lambda, \lambda, J, \phi) \qquad (B.1) $$

be a perturbed Hamiltonian system, and let $K(\Lambda, J, \phi)$ be a first integral of $H(\Lambda, \lambda, J, \phi)$ (such that $\{K, H\} = 0$, where $\{\cdot, \cdot\}$ is the usual Poisson bracket), and independent of $\lambda$. In addition, let

$$ H' = e^{L_s} H \qquad (B.2) $$

be a formal averaging of $H$ with respect to $\lambda$. If $S(\Lambda, \lambda, J, \phi)$ is a generator defined as below (Eqs. (B.4), (B.6), (B.7)), such that $H'$ is independent of $\lambda$. Then, $K$ is an integral of $H'$, i.e.

$$ \{K, H'\} = 0. \qquad (B.3) $$

The generator $S = \varepsilon S_1 + \varepsilon^2 S_2 + \varepsilon^3 S_3 + \dots$ is obtained formally through the identification order by order

$$ \begin{aligned} H'_0 &= H_0 \\ H'_1 &= H_1 + \{H_0, S_1\} \\ H'_2 &= \{H_0, S_2\} + \frac{1}{2}\{\{H_0, S_1\}, S_1\} + \{H_1, S_1\} \end{aligned} \qquad (B.4) $$

...

For any function $G(\Lambda, \lambda, J, \phi)$, let

$$ \langle G \rangle_{\lambda} = \frac{1}{(2\pi)^m} \int G d^m \lambda \qquad (B.5) $$

be the average of $G$ over all the angles $\lambda_k$. For each $n \ge 1$, the $S_n$ is obtained through the resolution of an equation of the form

$$ H'_n = \{H_0, S_n\} + \mathcal{R}_n, \qquad (B.6) $$

where $\mathcal{R}_n$ belongs to $\mathcal{L}(H_0, H_1, S_1, \dots, S_{n-1})$, the Lie algebra generated by $(H_0, H_1, S_1, \dots, S_{n-1})$. $H'_n$ will be the averaged part of $\mathcal{R}_n$, $\langle \mathcal{R}_n \rangle_{\lambda}$, and $S_n$ is obtained by solving the homological equation

$$ \{H_0, S_n\} = \langle \mathcal{R}_n \rangle_{\lambda} - \mathcal{R}_n. \qquad (B.7) $$

We show by recurrence that $\{K, S_n\} = 0$ for all $n \ge 1$. First, we note that as $\{K, H_0\} = 0$, we also have $\{K, H_1\} = 0$. As $K$ is independent of $\lambda_k$, we also have for all $G$

$$ \langle \{K, G\} \rangle_{\lambda} = \{K, \langle G \rangle_{\lambda}\}. \qquad (B.8) $$

This can be seen by formal expansion of $G$ in a Fourier series $G = \sum g_k e^{i(k, \lambda)}$. We have thus $\langle G \rangle_{\lambda} = g_0$. Let us assume now that $\{K, S_k\} = 0$ for all $k \le n$. As $\mathcal{R}_{n+1} \in \mathcal{L}(H_0, H_1, S_1, \dots, S_n)$, we also have $\{K, \mathcal{R}_{n+1}\} = 0$, and from Eq. (B.8) $\{K, \langle \mathcal{R}_{n+1} \rangle_{\lambda}\} = 0$ and thus

$$ \{K, \{H_0, S_{n+1}\}\} = 0. \qquad (B.9) $$

Using the Jacobi identity, as $\{F, H_0\} = 0$, we have

$$ \{H_0, \{K, S_{n+1}\}\} = 0. \qquad (B.10) $$

The solution of the homological Eq. (B.7) is unique up to a term independent of $\lambda$. But as $\langle \{K, S_{n+1}\} \rangle_{\lambda} = 0$, then the only possible solution for Eq. (B.10) is

$$ \{K, S_{n+1}\} = 0. \qquad (B.11) $$

In the same way, as $H'_1 = \langle H_1 \rangle_{\lambda}$, we have $\{K, H'_1\} = \langle \{K, H_1\} \rangle_{\lambda} = 0$. Thus $\{K, \{H_0, S_1\}\} = 0$, by the Jacobi identity, $\{H_0, \{K, S_1\}\} = 0$, and as previously, $\{K, S_1\} = 0$. Our recurrence is thus complete and it follows immediately that $\{K, H'\} = 0$.
# Appendix C: Special values of $C_c(\alpha, \gamma)$

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

## C.1. Asymptotic value of $C_c(\alpha, \gamma)$ for $\gamma \to 0$

We have shown that $e_c(\alpha, \gamma)$ is a differentiable function of $\gamma$, which is monotonic (41) and bounded ($e_c(\alpha, \gamma) \in [0, 1]$). The limit $e_c(\alpha, 0) = \lim_{\gamma \to 0} e_c(\alpha, \gamma)$ exists for all $\alpha \in ]0, 1]$. If $e_c(\alpha, \gamma)$ 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 $\alpha(1 + e)$:

$$ K(e, \alpha, \gamma) = \alpha(\gamma^2 - \alpha)e^3 - (2 - \alpha)(\gamma^2 - \alpha)e^2 - (1 - \alpha^2)e + (1 - \alpha)^2 = 0 \quad (C.1) $$

As $e_c(\alpha, \gamma)$, is a continuous function of $\gamma$, when $\gamma \to 0$ the limits $e_{0c}(\alpha)$ will satisfy the limit equation

$$ K_0(e, \alpha) = (1 - \alpha - \alpha e)^2 (1 - e) = 0 \quad (C.2) $$

with solutions $e_0 = 1/\alpha - 1$ and $e_1 = 1$. Depending on $\alpha$, several cases are treated:

### $\alpha < 1/2$: we have then $e_0 > 1$, and the only possibility for $e_c(\alpha, 0)$ is $e_1 = 1$; as it is the only root of Eq. (C.2) which belongs to $[0, 1]$, we have thus

$$ \lim_{\gamma \to 0} e_c(\alpha, \gamma) = 1. \quad (C.3) $$

We have then

$$ \begin{aligned} \lim_{\gamma \to 0} e'_{c}(\alpha, \gamma) &= 1 - 2\alpha; \\ \lim_{\gamma \to 0} C_c(\alpha, \gamma) &= 1 - 2\sqrt{\alpha(1 - \alpha)}. \end{aligned} \quad (C.4) $$

In order to study the behaviour of $e_c(\alpha, \gamma)$ in the vicinity of $\gamma = 0$, we can differentiate $K(e_c(\alpha, \gamma), \gamma) = 0$ twice, which gives

$$ \frac{de_c}{d\gamma}(\alpha, 0) = 0; \quad \frac{d^2 e_c}{d\gamma^2}(\alpha, 0) = \frac{4(\alpha - 1)}{(1 - 2\alpha)^2} < 0, \quad (C.5) $$

thus $e_c(\alpha, \gamma)$ decreases with respect to $\gamma$ at $\gamma = 0$. The second order development of $C_c$ gives

$$ C_c(\alpha, \gamma) = 1 - 2\sqrt{\alpha(1 - \alpha)} + \gamma\sqrt{\alpha} - \gamma^2\sqrt{\frac{1}{\alpha} - 1} + O(\gamma^3). \quad (C.6) $$

### $\alpha > 1/2$: in this case, $e_0 < 1$. As $e_c(\alpha, \gamma) \in ]0, e_0[$, we have $e_c(\alpha, 0) \in [0, e_0]$, which gives as the unique possibility

$$ \lim_{\gamma \to 0} e_c(\alpha, \gamma) = e_0 \quad (C.7) $$

with

$$ \lim_{\gamma \to 0} e'_{c}(\alpha, \gamma) = 0; \quad \lim_{\gamma \to 0} C_c(\alpha, \gamma) = 0. \quad (C.8) $$

By setting $\gamma = 0$ in Eq. (41), we also obtain

$$ \frac{de_c}{d\gamma}(\alpha, 0) = -\frac{e_0}{\alpha^{3/2}\sqrt{1 - e_0^2}} < 0. \quad (C.9) $$

The development of $C_c$ gives

$$ C_c(\alpha, \gamma) = \gamma \left( \sqrt{\alpha} - \sqrt{2 - \frac{1}{\alpha}} \right) - \frac{\gamma^2}{2\alpha - 1} \frac{(1 - \alpha)^2}{\alpha} + O(\gamma^3). \quad (C.10) $$

$\alpha = 1/2$: in this case, $e_0 = 1$, and the only solution is

$$ \lim_{\gamma \to 0} e_c(\alpha, \gamma) = 1 \quad (C.11) $$

and

$$ \lim_{\gamma \to 0} e'_{c}(\alpha, \gamma) = 0; \quad \lim_{\gamma \to 0} C_c(\alpha, \gamma) = 0. \quad (C.12) $$

Moreover, equation (C.1) becomes $\gamma^2 2e^2(3 - e) = (1 - e)^3$. We obtain thus the asymptotic value for $e_c(1/2, \gamma)$ when $\gamma \to 0$ as

$$ 1 - e_c \left( \frac{1}{2}, \gamma \right) \sim (4\gamma)^{2/3}. \quad (C.13) $$

For $\alpha = 1/2$, the development of $C_c$ in $\gamma$ contains non-polynomial terms in $\gamma$ giving

$$ C_c(1/2, \gamma) = \frac{\gamma}{\sqrt{2}} - 2^{-1/3}\gamma^{4/3} - 2^{-4/3}\gamma^2 + O(\gamma^{8/3}). \quad (C.14) $$

## C.2. Asymptotic value of $C_c(\alpha, \gamma)$ for $\gamma \to +\infty$

This case is more simple. If $e_c(\alpha, \gamma)$ is a solution of Eq. (35), then it will also be a solution of Eq. (C.1), and thus of

$$ \frac{K(e, \alpha, \gamma)}{\gamma^2} = 0. \quad (C.15) $$

As $e_c(\alpha, \gamma)$ is monotonic and bounded, it has a limit when $\gamma \to +\infty$, which will verify the limit equation (C.15), when $\gamma \to +\infty$, that is

$$ K_{\infty}(\alpha, e) = e^2(2 - \alpha - \alpha e) = 0. \quad (C.16) $$

As $0 < \alpha < 1$, the only solution is $e = 0$, and thus

$$ \lim_{\gamma \to +\infty} e_c(\alpha, \gamma) = 0 \quad (C.17) $$

and

$$ \lim_{\gamma \to +\infty} e'_{c}(\alpha, \gamma) = 1 - \alpha; \quad \lim_{\gamma \to +\infty} C_c(\alpha, \gamma) = 1 - \sqrt{\alpha(2 - \alpha)} \quad (C.18) $$

and more precisely

$$ C(\alpha, \gamma) = 1 - \sqrt{\alpha(2 - \alpha)} - \frac{1}{2\gamma} \frac{\sqrt{\alpha(1 - \alpha)^2}}{2 - \alpha} + \frac{1}{\gamma^2} \frac{\sqrt{\alpha(2 - \alpha)(1 - \alpha)^2}}{2(2 - \alpha)^3} + O\left(\frac{1}{\gamma^3}\right). \quad (C.19) $$

## C.3. Study of $C_c(\alpha, \gamma)$ for $\gamma = 1$ and $\gamma = \sqrt{\alpha}$

For $\gamma = 1$ or $\gamma = \sqrt{\alpha}$, we can also obtain simple expressions for $C_c(\alpha, \gamma)$. Indeed, If $\gamma = 1$, $K(e, \alpha, \gamma)$ factorises in $(1 - \alpha)(1 + e)(\alpha e^2 - 2e + 1 - \alpha)$ and we have a single solution for $e_c$ in the interval $[0, 1]$,

$$ e_c(\alpha, 1) = \frac{1 - \sqrt{1 - \alpha + \alpha^2}}{\alpha}; \quad (C.20) $$

$$ e'_{c}(\alpha, 1) = \sqrt{1 - \alpha + \alpha^2} - \alpha; $$
and

$$C_c(\alpha, 1) = \sqrt{\alpha} - \frac{\sqrt{\alpha - 2 + 2\sqrt{1 - \alpha + \alpha^2}}}{\sqrt{\alpha}} \quad (C.21) \\ +1 - \sqrt{\alpha} \sqrt{1 - 2\alpha + 2\sqrt{1 - \alpha + \alpha^2}}.$$

For $\gamma = \sqrt{\alpha}$, the cubic Eq. (C.1) reduces to

$$K(\alpha, \sqrt{\alpha}) = -(1 - \alpha^2)e + (1 - \alpha)^2 = 0 \quad (C.22)$$

with the single solution

$$e_c(\alpha, \sqrt{\alpha}) = e'_c(\alpha, \sqrt{\alpha}) = \frac{1 - \alpha}{1 + \alpha} \quad (C.23)$$

and

$$C_c(\alpha, \sqrt{\alpha}) = (1 - \sqrt{\alpha})^2. \quad (C.24)$$

### C.4. $C_c(\alpha, \gamma)$ for $\alpha \to 1$

Let us denote $\eta = 1 - \alpha$. Equation (C.1) can be developed in $\eta$,

$$\begin{aligned} K(e, \alpha, \gamma) &= (\gamma^2 - 1)e^2(e - 1) \\ &= \eta e((\gamma^2 - 1)(e - e^2) - 2) \\ &= \eta^2(1 - e - e^2 - e^3) = 0. \end{aligned} \quad (C.25)$$

The zeroth and first orders of Eq. (C.25) imply that $e_c$ must go to zero; moreover, it scales with $\eta$. We write $e_c(\eta, \gamma) = \kappa(\gamma)\eta + o(\eta)$. We inject this expression in Eq. (C.25) and keep the second order in $\eta$

$$(\gamma^2 - 1)\kappa^2 + 2\kappa - 1 = 0. \quad (C.26)$$

We keep the solution that is positive and continuous in $\gamma$ and we have

$$e_c(\eta, \gamma) = \frac{\eta}{\gamma + 1} + o(\eta). \quad (C.27)$$

If we now compute $C_c$ developed for $\alpha \to 1$ we have

$$\begin{aligned} C_c(1 - \eta, \gamma) &= k(\gamma)\eta^2 + o(\eta^2) \\ &= \frac{1}{2} \frac{\gamma}{\gamma + 1} \eta^2 + o(\eta^2). \end{aligned} \quad (C.28)$$