AMD-stability in presence of first order Mean Motion Resonances

TL;DR

This paper extends the AMD-stability criterion to include first order mean-motion resonances, providing a unified stability criterion that accounts for resonance overlap and short-term chaos in planetary systems.

Contribution

It introduces a new overlap criterion for first order MMR and refines AMD-stability to incorporate chaos from resonance overlap, unifying previous stability criteria.

Findings

Derived a new resonance overlap criterion for first order MMR.

Unified existing stability criteria into a comprehensive framework.

Applied the improved AMD-stability to real multi-planet systems.

Abstract

The AMD-stability criterion allows to discriminate between a-priori stable planetary systems and systems for which the stability is not granted and needs further investigations. AMD-stability is based on the conservation of the Angular Momentum Deficit (AMD) in the averaged system at all orders of averaging. While the AMD criterion is rigorous, the conservation of the AMD is only granted in absence of mean-motion resonances (MMR). Here we extend the AMD-stability criterion to take into account mean-motion resonances, and more specifically the overlap of first order MMR. If the MMR islands overlap, the system will experience generalized chaos leading to instability. The Hamiltonian of two massive planets on coplanar quasi-circular orbits can be reduced to an integrable one degree of freedom problem for period ratios close to a first order MMR. We use the reduced Hamiltonian to derive a new overlap criterion for first order MMR. This stability criterion unifies the previous criteria proposed in the literature and admits the criteria obtained for initially circular and eccentric orbits as limit cases. We then improve the definition of AMD-stability to take into account the short term chaos generated by MMR overlap. We analyze the outcome of this improved definition of AMD-stability on selected multi-planet systems from the Extrasolar Planets Encyclopeadia.

Tables

Table 1. Table 1: Summary of the diverse notations of AMD used in this paper.

Notation	Description	Equation
$C$	Total AMD of the system
$C_{\min }$	Minimal AMD to enter a resonance island	(62)
$c_{\min }$	Normalized minimal AMD	(65)
$𝒞$	Relative AMD	(91)
$C_c^{\left(0\right)}$	Critical AMD deduced from the collision condition	(Laskar & Petit, 2017)
$C_c^{\left(1\right)}$	Critical AMD deduced from the MMR overlap	(97)
$C_c$	Complete critical AMD	(101)

Table 2. Table E.1: AMD-stability coefficients computed for the systems affected by the MMR overlap criterion

Planet	Period (d)	Mass (${ℳ}_E^N$)	Eccentricity	$\sqrt{⟨e^2⟩}$	$\beta$	${\beta }^{\left(\mathrm{MMR}\right)}$
HD 128311	Mass: 0.84 $M_{\odot }^N$
b	454.2	$463.14$	$0.345$	$0.352$	$0.312$
c	923.8	$1032.46$	$0.230$	$0.244$	$3.200$	$27.931$
HD 200964	Mass: 1.44 $M_{\odot }^N$
b	613.8	$587.98$	$0.040$	$0.067$	$0.024$
c	825	$284.46$	$0.181$	$0.184$	$3.872$	$+\infty$
HD 204313	Mass: 1.045 $M_{\odot }^N$
c	34.905	$17.58$	$0.155$	$0.184$	$16.664$
b	2024.1	$1360.31$	$0.095$	$0.095$	$0.110$	$0.110$
d	2831.6	$533.95$	$0.280$	$0.308$	$8.032$	$+\infty$
HD 33844	Mass: 1.75 $M_{\odot }^N$
b	551.4	$622.94$	$0.150$	$0.180$	$0.084$
c	916	$556.20$	$0.130$	$0.189$	$2.939$	$22.676$
HD 45364	Mass: 0.82 $M_{\odot }^N$
b	226.93	$59.50$	$0.168$	$0.171$	$0.070$
c	342.85	$209.10$	$0.097$	$0.099$	$1.975$	$13.700$
HD 47366	Mass: 1.81 $M_{\odot }^N$
b	363.3	$556.20$	$0.089$	$0.138$	$0.146$
c	684.7	$591.16$	$0.278$	$0.292$	$2.896$	$2.896$
HD 5319	Mass: 1.56 $M_{\odot }^N$
b	675	$616.59$	$0.120$	$0.162$	$0.053$
c	886	$365.50$	$0.150$	$0.171$	$8.659$	$+\infty$
HD 73526	Mass: 1.08 $M_{\odot }^N$
b	188.9	$715.11$	$0.290$	$0.293$	$0.200$
c	379.1	$715.11$	$0.280$	$0.289$	$3.391$	$9.922$

Equations

$$\hat{\mathcal{H}}=$$

$$\frac{1}{2}\frac{\Vert {\tilde{\mathbf{u}}}_1+{\tilde{\mathbf{u}}}_2{\Vert }^2}{m_0}-\mathcal{G}\frac{m_1{m}_2}{{\Delta }_{12}}$$

$$\hat{\mathcal{K}}$$

$$\epsilon {\hat{\mathcal{H}}}_1$$

$$\epsilon =\frac{m_1+m_2}{m_0}.$$

$$\text{@vec}\hat{G}=\sum _{i=1}^2{\mathbf{u}}_i\wedge {\tilde{\mathbf{u}}}_i$$

$$\hat{\mathcal{H}}$$

$$=-\sum _{i=1}^2\frac{{\mu }^2{m}_i^3}{2{\hat{\Lambda }}_i^2}+\epsilon \underset{k\in {\mathbb{Z}}^2}{\sum _{l,\bar{l}\in {\mathbb{N}}^2}}{C}_{l,\bar{l},k}(\hat{\Lambda })\prod _{i=1}^2{\hat{x}}_i^{l_i}{\bar{\hat{x}}}_i^{{\bar{l}}_i}{\mathbf{e}}^{\mathbf{i}{k}_i{\lambda }_i},$$

$${\hat{\Lambda }}_i=m_i\sqrt{\mu a_i}$$

$${\hat{G}}_i={\hat{\Lambda }}_i\sqrt{1-e_i^2}$$

$${\hat{C}}_i={\hat{\Lambda }}_i-{\hat{G}}_i$$

$${\hat{x}}_i=\sqrt{{\hat{C}}_i}{\mathbf{e}}^{-\mathbf{i}{\varpi }_i}$$

$${\lambda }_i=M_i+{\varpi }_i.$$

$$\sum _{i=1}^2{k}_i-l_i+{\bar{l}}_i=0.$$

$$\left(\begin{array}[]{c}\sigma\\ M_{2}\end{array}\right)=\left(\begin{array}[]{cc}-p&p+1\\ 0&1\end{array}\right)\left(\begin{array}[]{c}M_{1}\\ M_{2}\end{array}\right).$$

$$\left(\begin{array}[]{c}\hat{\Gamma}_{1}\\ \hat{\Gamma}\end{array}\right)=\left(\begin{array}[]{cc}-{{1}\over{p}}&0\\ {{p+1}\over{p}}&1\end{array}\right)\left(\begin{array}[]{c}\hat{\Lambda}_{1}\\ \hat{\Lambda}_{2}\end{array}\right)=\left(\begin{array}[]{c}{{-1}\over{p}}\hat{\Lambda}_{1}\\ {{p+1}\over{p}}\hat{\Lambda}_{1}+\hat{\Lambda}_{2}\end{array}\right).$$

$$\hat{\Gamma }=\frac{p+1}{p}{\hat{\Lambda }}_1+{\hat{\Lambda }}_2$$

$${\hat{\mathcal{H}}}_{\mathrm{av}}=$$

$$\epsilon \underset{j\in \mathbb{Z}}{\sum _{l,\bar{l}\in {\mathbb{N}}^2}}{C}_{l,\bar{l},j}(\hat{\Lambda }){\hat{x}}_1^{l_1}{\bar{\hat{x}}}_1^{{\bar{l}}_1}{\hat{x}}_2^{l_2}{\bar{\hat{x}}}_2^{{\bar{l}}_2}{\mathbf{e}}^{\mathbf{i}j((p+1){\lambda }_2-p{\lambda }_1)},$$

$$\left(\begin{array}[]{c}\theta_{\Gamma}\\ \theta_{G}\\ \sigma_{1}\\ \sigma_{2}\end{array}\right)=\left(\begin{array}[]{cccc}p&-p&0&0\\ -p&p+1&0&0\\ -p&p+1&1&0\\ -p&p+1&0&1\end{array}\right)\cdot\left(\begin{array}[]{c}\lambda_{1}\\ \lambda_{2}\\ -\varpi_{1}\\ -\varpi_{2}\end{array}\right).$$

$$\left(\begin{array}[]{c}\hat{\Gamma}\\ \hat{G}\\ \hat{C}_{1}\\ \hat{C}_{2}\end{array}\right)=\left(\begin{array}[]{cccc}{{p+1}\over{p}}&1&0&0\\ 1&1&-1&-1\\ 0&0&1&0\\ 0&0&0&1\end{array}\right)\cdot\left(\begin{array}[]{c}\hat{\Lambda}_{1}\\ \hat{\Lambda}_{2}\\ \hat{C}_{1}\\ \hat{C}_{2}\end{array}\right).$$

$$\prod _{i=1}^2{\hat{x}}_i^{l_i}{\bar{\hat{x}}}_i^{{\bar{l}}_i}{\mathbf{e}}^{\mathbf{i}j{\theta }_G}$$

$$=\prod _{i=1}^2{\hat{\mathfrak{X}}}_i^{l_i}{\overline{\hat{\mathfrak{X}}}}_i^{{\bar{l}}_i};$$

$${\hat{\mathcal{H}}}_{\mathrm{av}}=-\sum _{i=1}^2\frac{{\mu }^2{m}_i^3}{2{\hat{\Lambda }}_i^2}+\epsilon \underset{j\in \mathbb{Z}}{\sum _{l,\bar{l}\in {\mathbb{N}}^2}}{C}_{l,\bar{l},j}(\hat{\Lambda })\prod _{i=1}^2{\hat{\mathfrak{X}}}_i^{l_i}{\overline{\hat{\mathfrak{X}}}}_i^{{\bar{l}}_i}.$$

$${\hat{\Lambda }}_1$$

$${\hat{\Lambda }}_2$$

$$\alpha =\frac{a_1}{a_2},$$

$$\epsilon {\hat{\mathcal{H}}}_1={\hat{R}}_1({\hat{\mathfrak{X}}}_1+{\overline{\hat{\mathfrak{X}}}}_1)+{\hat{R}}_2({\hat{\mathfrak{X}}}_2+{\overline{\hat{\mathfrak{X}}}}_2),$$

$${\hat{R}}_1$$

$$\mathrm{and}$$

$${\hat{R}}_2$$

$$r_1(\alpha )$$

$$\mathrm{and}$$

$$r_2(\alpha )$$

$$+\frac{1}{2}{b}_{1/2}^{(p)}(\alpha ).$$

$$b_s^{(k)}(\alpha )=\frac{1}{\pi }{\int }_{-\pi }^{\pi }\frac{\cos (k\varphi )}{{\left(1-2\alpha \cos \varphi +{\alpha }^2\right)}^s}\mathrm{d}\varphi$$

$${\mathcal{H}}_{1,i}=\frac{{\mu }^2{m}_1^2{m}_2^2}{2m_0{\hat{\Lambda }}_1{\hat{\Lambda }}_2}\sqrt{\frac{2}{{\hat{\Lambda }}_2}}({\hat{\mathfrak{X}}}_2+{\overline{\hat{\mathfrak{X}}}}_2).$$

$${\alpha }_0={\left(\frac{p}{p+1}\right)}^{2/3},$$

Full text

1 Introduction

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

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

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

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

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

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

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

2 The resonant Hamiltonian

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

2.1 Averaged Hamiltonian in the vicinity of a resonance

Let us consider two planets of masses $m_{1}$ and $m_{2}$ orbiting a star of mass $m_{0}$ in the plane. We denote the positions of the planets, ${\bf u}_{i}$, and the associated canonical momenta in the heliocentric frame, ${\bf\tilde{u}}_{i}$. The Hamiltonian of the system is (Laskar & Robutel, 1995)

[TABLE]

where $\Delta_{12}=\|{\bf u}_{1}-{\bf u}_{2}\|$, and $\mathcal{G}$ is the constant of gravitation. $\hat{\mathcal{H}}$ can be decomposed into a Keplerian part $\hat{\mathcal{K}}$ describing the motion of the planets if they had no masses and a perturbation part $\varepsilon\hat{\mathcal{H}}_{1}$ due to the influence of massive planets,

[TABLE]

The small parameter $\varepsilon$ is defined as the ratio of the planet masses over the star mass

[TABLE]

Let us denote the angular momentum,

[TABLE]

which is simply the sum of the two planets Keplerian angular momentum. $\@vec{\hat{G}}$ is a first integral of the system.

Following (Laskar, 1991), we express the Hamiltonian in terms of the Poincaré coordinates

[TABLE]

where $\mu=\mathcal{G}m_{0}$ and for $i=1,2$,

[TABLE]

Here, $M_{i}$ corresponds to the mean anomaly, $\varpi_{i}$ to the longitude of the pericenter, $a_{i}$ to the semi-major axis and $e_{i}$ to the eccentricity of the Keplerian orbit of the planet $i$. $\hat{G}_{i}$ is the Keplerian angular momentum of planet $i$. We use the set of symplectic coordinates of the problem $(\hat{\Lambda}_{i},\lambda_{i},\hat{C}_{i},-\varpi)$, or the canonically associated variables $(\hat{\Lambda}_{i},\lambda_{i},\hat{x}_{i},-\mathbf{i}\bar{\hat{x}}_{i})$. The coefficients $C_{l,\bar{l},k}$ depend on $\hat{\Lambda}$ and the masses of the bodies. They are linear combinations of Laplace coefficients (Laskar & Robutel, 1995). As a consequence of angular momentum conservation, the d’Alembert rule gives a relation on the indices of the non-zero $C_{l,\bar{l},k}$ coefficients

[TABLE]

We study here a system with periods close to the first-order MMR $p:p+1$ with $p\in\mathbb{N}^{*}$. For periods close to this configuration, we have ${-pn_{1}+(p+1)n_{2}\simeq 0}$, where $n_{i}=\mu^{2}m_{i}^{3}/\hat{\Lambda}_{i}^{3}$ is the Keplerian mean motion of the planet $i$.

2.1.1 Averaging over non-resonant mean-motions

Due to the $p:p+1$ resonance, we cannot average on both mean anomalies independently. Therefore, there is no conservation of $\hat{\Lambda}_{i}$ as in the secular problem. However, the partial averaging over one of the mean anomaly gives another first integral. Following (Delisle et al., 2012), we consider the equivalent set of coordinates $(\hat{\Lambda}_{i},M_{i},\hat{G}_{i},\varpi_{i})$, and make the following change of angles

[TABLE]

The actions associated to these angles are

[TABLE]

We can now average the Hamiltonian over $M_{2}$ using a change of variables close to the identity given by the Lie series method. Up to terms of orders $\varepsilon^{2}$, we can kill all the terms with indices not of the form $C_{l,\bar{l},-jp,j(p+1)}$. In order to keep the notations light, we do not change the name of the variables after the averaging. We also designate the remaining coefficients $C_{l,\bar{l},-jp,j(p+1)}$ by the lighter expression $C_{l,\bar{l},j}$. Since $M_{2}$ does not appear explicitly in the remaining terms,

[TABLE]

is a first integral of the averaged Hamiltonian. The parameter $p\hat{\Gamma}$ is often designed as the spacing parameter (Michtchenko et al., 2008) and has been used extensively in the study of the first-order MMR dynamics.

Expressed with the variables $(\hat{\Lambda},\lambda,\hat{x},\bar{\hat{x}})$, the Hamiltonian can be written

[TABLE]

where we dropped the terms of order $\varepsilon^{2}$.

2.1.2 Poincare-like complex coordinates

Delisle et al. (2012) used a change of the angular coordinates in order to remove the exponential in the second term of eq. (11) and use $\hat{G}$ and $\hat{\Gamma}$ as actions. The new set of angles $(\theta_{\Gamma},\theta_{G},\sigma_{1},\sigma_{2})$ is defined as

[TABLE]

The conjugated actions are

[TABLE]

We define $\hat{\mathfrak{\scriptstyle X}}_{i}=\sqrt{\hat{C}_{i}}\mathbf{e}^{\mathbf{i}\sigma_{i}}$, the complex coordinates associated to $(\hat{C}_{i},\sigma_{i})$. Since we have $\hat{\mathfrak{\scriptstyle X}}_{i}=\hat{x}_{i}\mathbf{e}^{\mathbf{i}\theta_{G}}$, the terms of the perturbation in (11) can be written

[TABLE]

the last equality resulting from the d’Alembert rule (7). $\hat{\Gamma}$ and $\hat{G}$ are conserved and the averaged Hamiltonian no longer depends on the angles $\theta_{\Gamma}$ and $\theta_{G}$

[TABLE]

$\hat{\Lambda}_{1}$ and $\hat{\Lambda}_{2}$ can be expressed as functions of the new variables and we have

[TABLE]

where $\hat{C}=\hat{C}_{1}+\hat{C}_{2}$ is the total AMD of the system. Up to the value of the first integrals $\hat{\Gamma}$ and $\hat{G}$, the system now has two effective degrees of freedom.

2.2 Computation of the perturbation coefficients

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

We follow the method described in Laskar (1991) and Laskar & Robutel (1995) to determine the expression of the perturbation $\hat{\mathcal{H}}_{1}$. The details of the computation are given in Appendix A. Since we compute an expression at first order in eccentricities and $\varepsilon$, the semi major axis and in particular their ratio,

[TABLE]

are evaluated at the resonance. At the first order, the perturbation term $\hat{\mathcal{H}}_{1}$ has for expression

[TABLE]

where

[TABLE]

with $\gamma=m_{1}/m_{2}$,

[TABLE]

In the two previous expressions, $b^{(k)}_{s}(\alpha)$ are the Laplace coefficients that can be expressed as

[TABLE]

for $k>0$. For $k=0$, a $1/2$ factor has to be added in the second-hand member of (25).

For $p=1$, it should be noted that a contribution from the kinetic part should be added (Appendix A and Delisle et al., 2012)

[TABLE]

Using the expression of $\alpha$ at the resonance ${p:p+1}$,

[TABLE]

we can give the asymptotic development of the coefficients $r_{1}$ and $r_{2}$ for $p\rightarrow+\infty$ (see Appendix A.1). The equivalent is

[TABLE]

where $K_{\nu}(x)$ is the modified Bessel function of the second kind. We note $r$ the numerical factor of the equivalent (28), we have

[TABLE]

For the resonant coefficients $r_{1}$ and $r_{2}$, Deck et al. (2013) used the expressions $f_{p+1,27}(\alpha)$ and $f_{p+1,31}(\alpha)$ given in (Murray & Dermott, 1999, pp. 539-556). The expressions (23) and (24) are similar to $f_{p+1,27}(\alpha)$ and $f_{p+1,31}(\alpha)$ up to algebraic transformations using the relations between Laplace coefficients (Laskar & Robutel, 1995). In their computations, Deck et al. used a numerical fit of the coefficients for $p=2$ to 150 and obtained

[TABLE]

We obtain the same numerical factor $r$ through the analytical development of the functions $r_{1}$ and $r_{2}$.

2.3 Renormalization

So far, the Hamiltonian has two degrees of freedom $(\hat{\mathfrak{\scriptstyle X}}_{1},\bar{\hat{\mathfrak{\scriptstyle X}}}_{1},\hat{\mathfrak{\scriptstyle X}}_{2},\bar{\hat{\mathfrak{\scriptstyle X}}}_{2})$ and depends on two parameters $\hat{G}$ and $\hat{\Gamma}$. As shown in (Delisle et al., 2012), the constant $\hat{\Gamma}$ can be used to scale the actions, the Hamiltonian and the time without modifying the dynamics. We define

[TABLE]

With this change of variables, the new Hamiltonian no longer depends on $\hat{\Gamma}$.

The shape of the phase space is now only dependent on the first integral $G$. However, $G$ does not vanish for the configuration around which the Hamiltonian is developed: the case of two resonant planets on circular orbits. To be able to develop the Keplerian part in power of the system’s parameter, we define $\Delta G=G_{0}-G$, the difference in angular momentum between the circular resonant system and the actual configuration. We have

[TABLE]

where $\Lambda_{1,0}$ and $\Lambda_{2,0}$ are the value of $\Lambda_{1}$ and $\Lambda_{2}$ at resonance. By definition, we have

[TABLE]

Moreover, we can express $\Lambda_{1,0}$ as a function of the ratios $\alpha_{0}$ and $\gamma$,

[TABLE]

Similarly, $\Lambda_{2,0}$ can be expressed as

[TABLE]

Since $G_{0}$ is constant, $\Delta G$ is also a first integral of $\mathcal{H}$. From now on, we consider $\Delta G$ as a parameter of the two-degrees-of-freedom $(\mathfrak{\scriptstyle X}_{1},\mathfrak{\scriptstyle X}_{2})$ Hamiltonian $\mathcal{H}$. The Keplerian part depends on the coordinates $\mathfrak{\scriptstyle X}_{i}$ through the dependence of $\Lambda_{i}$ in $C$.

$\Lambda_{1}$ and $\Lambda_{2}$ can be expressed as functions of the Hamiltonian coordinates and their value at the resonance,

[TABLE]

2.4 Integrable Hamiltonian

The system can be made integrable by a rotation of the coordinates $\mathfrak{\scriptstyle X}_{i}$ (Sessin & Ferraz-Mello, 1984; Henrard et al., 1986; Delisle et al., 2014). We introduce $R$ and $\phi$ such that

[TABLE]

We have $R^{2}=R_{1}^{2}+R_{2}^{2}$ and $\tan(\phi)=R_{2}/R_{1}$. If we note $\mathcal{R}_{\phi}$ the rotation of angle $\phi$ we define $y$ such that $\mathfrak{\scriptstyle X}=\mathcal{R}_{\phi}y$. We still have $C=\sum y_{i}\bar{y}_{i}$ so the only change in the Hamiltonian is the perturbation term

[TABLE]

where $(I,\theta)$ are the action-angle coordinates associated to $y$. With these coordinates, $I_{2}$ is a first integral. $R$ has for expression

[TABLE]

We now develop the Keplerian part around the circular resonant configuration in series of $(C-\Delta G)$ thanks to the relations (35). We develop the Keplerian part to the second order in $(C-\Delta G)$ since the first order vanishes (see Appendix B). The computation of the second-order coefficient gives

[TABLE]

We drop the constant part of the Hamiltonian and obtain the following expression

[TABLE]

We again change the time scale by dividing the Hamiltonian by $-\mathcal{K}_{2}$ and multiplying the time by this factor. We define

[TABLE]

and after simplification,

[TABLE]

where $r$ was defined in (29) and ${f(p)=1+\mathrm{O}(p^{-1})}$ is a function of $p$ and $\gamma$

[TABLE]

At this point the Hamiltonian can be written

[TABLE]

and has almost its final form. We divide the actions and the time by $\chi^{2/3}$ and the Hamiltonian by $\chi^{4/3}$ and we obtain

[TABLE]

where

[TABLE]

2.5 Andoyer Hamiltonian

We now perform a polar to Cartesian change of coordinates with

[TABLE]

We change the sign of $X$ in order to have the same orientation as (Deck et al., 2013). Doing so, the Hamiltonian becomes

[TABLE]

We recognize the second fundamental model of resonance (Henrard & Lemaitre, 1983). This Hamiltonian is also called an Andoyer Hamiltonian (Ferraz-Mello, 2007). We show in Figure 1 the level curves of the Hamiltonian $\mathcal{H}_{A}$ for $\mathcal{I}_{0}=3$.

The fixed points of the Hamiltonian satisfy the equations

[TABLE]

which have for solutions $Y=0$ and the real roots of the cubic equation in $X$

[TABLE]

Equation (52) has three solutions (Deck et al., 2013) if its determinant $\Delta=32(\mathcal{I}_{0}^{3}-27/8)>0$, i.e. $\mathcal{I}_{0}>3/2$. In this case, we note these roots $X_{1}<X_{2}<X_{3}$. $X_{1}$ and $X_{2}$ are elliptic fixed points while $X_{3}$ is a hyperbolic one.

3 Overlap criterion

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

Wisdom (1980) first applied the resonance overlap criterion to the first-order MMR and found, in the case of the restricted three-body problem with a circular planet, that the overlap occurs for

[TABLE]

Through numerical simulations, (Duncan et al., 1989) confirmed Wisdom’s expression up to the numerical coefficient ($1-\alpha<1.5\ \varepsilon^{2/7}$). A similar criterion was then developed by Mustill & Wyatt (2012) for an eccentric planet, they found that for an eccentricity above $0.2\ \varepsilon^{3/7}$, the overlap region satisfies the criterion $1-\alpha<1.8(\varepsilon\ e)^{1/5}$. Deck et al. (2013) adapted those two criteria to the case of two massive planets, finding little difference up to the numerical coefficients. However, they treat two different situations; the case of orbits initially circular and the case of two eccentric orbits. As in (Mustill & Wyatt, 2012), the eccentric criterion proposed can be used for eccentricities verifying $e_{1}+e_{2}\gtrsim 1.33\ \varepsilon^{3/7}$. We show here that the two Deck’s criteria can be obtained as the limit cases of a general expression.

3.1 Width of the libration area

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

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

We note $X_{1}^{*},X_{2}^{*}$ the abscissas of the intersections between the separatrices and the $Y=0$ axis. Relations between $X_{1}^{*},X_{2}^{*}$, and $X_{3}$ can be derived (see Appendix C.1) and we obtain the expressions of $X^{*}_{1}$ and $X^{*}_{2}$ as functions of $X_{3}$ (Ferraz-Mello, 2007; Deck et al., 2013). We have

[TABLE]

The width of the libration zone $\delta X$ depends solely on the value of $X_{3}$,

[TABLE]

In order to study the overlap of resonance islands, we need the width of the resonance in terms of $\alpha$. Let us invert the previous change of variables in order to express the variation of $\alpha$ in terms of the variation of $X$. In this subsection, for any function $Q(X)$, we note

[TABLE]

The computation of $\delta\mathcal{I}$ (48) is straightforward from the computation of $\delta X$

[TABLE]

We then directly deduce $\delta I_{1}=\chi^{2/3}\delta\mathcal{I}$ from (47). Since $I_{2}$ and $\Delta G$ are first integrals, the variation of $\Lambda_{i}$ only depends on $\delta I_{1}$. And finally, since we have

[TABLE]

$\alpha$ can be developed to the first order in $(C-\Delta G)$ thanks to (35). This development gives

[TABLE]

The width of the resonance in terms of $\alpha$ is then directly related to $X_{3}$ through

[TABLE]

The computation of the width of resonance is thus reduced to the computation of the root $X_{3}$ as a function of the parameters. It should also be remarked that at the first order, the width of resonance does not depend on the mass ratio $\gamma$.

3.2 Minimal AMD of a resonance

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

For a configuration close to a given resonance $p:p+1$, the AMD can evolve toward higher values if the original value places the system in a configuration above the inner separatrix, eventually leading the planets to collision or chaotic motion in case of MMR overlap. On the other hand, if the initial AMD of the planets forces them to remain in the inner circulation region of the overlapped MMR islands, the system will remain stable in regards to this criterion. Since $C=I_{1}+I_{2}$, and $I_{2}$ is a first integral, we define the minimal AMD of a resonance111We summarize the notations of the various AMD expressions used in this paper in Table 1. ${C_{\mathrm{min}}}(p)$ as the minimal value of $I_{1}$ to enter the resonant island given $\Delta G-I_{2}$. Two cases must be discussed:

• •

The point $I_{1}=0$ is already in the libration zone and then ${C_{\mathrm{min}}}=0$,

• •

The point $I_{1}=0$ is in the inner circulation zone and then we have

[TABLE]

In the second case, we have an implicit expression of $X_{3}$ depending on ${C_{\mathrm{min}}}$

[TABLE]

where $\chi$ was defined in (41). In other words, there is a one-to-one correspondence between ${C_{\mathrm{min}}}$ (62) and the Hamiltonian parameter $\mathcal{I}_{0}$ for ${C_{\mathrm{min}}}>0$. The shape of the resonance island is completely described by ${C_{\mathrm{min}}}$.

We can also use the definition of ${C_{\mathrm{min}}}$ to give an expression depending on the system parameters

[TABLE]

where $c_{i}=\sqrt{2}\sqrt{1-\sqrt{1-e_{i}^{2}}}=|\mathbf{\mathfrak{X}}_{i}|$. We note

[TABLE]

the reduced minimal AMD. We can use the expression (64) to compute the quantity $\chi^{-1/3}\sqrt{2{C_{\mathrm{min}}}}$ appearing in equation (63)

[TABLE]

The function ${C_{\mathrm{min}}}(X_{3})$ (Eq. 62) is plotted in Figure 2 with the two approximations used by Deck et al. (2013) to obtain the width of the resonance. For ${C_{\mathrm{min}}}\gg\chi^{2/3}$ or ${C_{\mathrm{min}}}$ close to zero, the relation can be simplified and we obtain

[TABLE]

We can use the developments (67) and (68) in order to compute the width of the resonance in these two cases (see Appendix C). It should be noted as well that for ${C_{\mathrm{min}}}=0$, we have $X_{3}=2^{2/3}$.

3.3 Implicit overlap criterion

The overlap of MMR can be determined by finding the first integer $p$ such that the sum of the half-width of the resonances $p:p+1$ and $p+1:p+2$ is larger than the distance between the respective centers of these two resonances (Wisdom, 1980; Deck et al., 2013)

[TABLE]

where $\Delta\alpha$ is the distance between the two centers and $\delta\alpha_{k}$ corresponds to the width of the resonance $k:k+1$.

Up to terms of order $\varepsilon^{2/3}$, the center of the resonance island $p:p+1$ is located at the center of the resonance of the unperturbed Keplerian problem, ${\alpha_{0,p}=(p/(p+1))^{2/3}}$. We develop $\alpha_{0,p}$ for $p\gg 1$

[TABLE]

Therefore, we have at second order in $p$

[TABLE]

We can use the implicit expression (63) of $X_{3}$ as a function of $\sqrt{c_{\mathrm{min}}}$ (Eq. 65) in order to derive an overlap criterion independent of approximations on the value of ${C_{\mathrm{min}}}$. Equating the general width of resonance (61) with the distance between to adjacent centers (71) and isolating $X_{3}$ gives

[TABLE]

We can inject this expression of $X_{3}$ into (63), and using equation (66),

[TABLE]

Using the first order expression of $(p+1)$ as a function of $\alpha$,

[TABLE]

we obtain an implicit expression of the overlap criterion

[TABLE]

3.4 Overlap criterion for circular orbits

The implicit expression (75) can be used to find the criteria proposed by Deck et al. (2013) for circular and eccentric orbits. Let us first obtain the circular criterion by imposing $c_{\mathrm{min}}=0$ in equation (75)

[TABLE]

We can express $1-\alpha$ as a function of $\varepsilon$ and we obtain

[TABLE]

The exponent 2/7 was first proposed by Wisdom (1980) and the numerical factor $1.46$ is similar to the one found by Deck et al. (2013).

3.5 Overlap criterion for high-eccentricity orbits

For large eccentricity, Deck et al. (2013) proposes a criterion based on the development (67) of equation (63). This criterion is obtained from (75) by ignoring the second term of the right-hand side which leads to

[TABLE]

Isolating $1-\alpha$ gives

[TABLE]

This result is also similar to Deck’s one. For small $c_{\mathrm{min}}$, the criterion (79) is less restrictive than the criterion (77) obtained for circular orbits. The comparison of these two overlap criteria provides a minimal value of $c_{\mathrm{min}}$ for the validity of the eccentric criterion

[TABLE]

3.6 Overlap criterion for low-eccentricity orbits

For smaller eccentricities, we can develop the equation (75) for small $\sqrt{c_{\mathrm{min}}}$ and $\alpha$ close to ${\alpha_{\mathrm{cir}}=1-1.46\varepsilon^{2/7}}$, the critical semi major axis ratio for the circular overlap criterion (77). We have

[TABLE]

We develop the right-hand side at the first order in $(\alpha_{\mathrm{cir}}-\alpha)$ and evaluate the left-hand side for $\alpha=\alpha_{\mathrm{cir}}$ and after some simplifications obtain

[TABLE]

We inject the expression of $\alpha_{\mathrm{cir}}$ into this equation and obtain the following development of the overlap criterion for low eccentricity:

[TABLE]

This development remains valid for small enough $\sqrt{c_{\mathrm{min}}}$ if ${\alpha_{\mathrm{cir}}-\alpha\ll 1-\alpha_{\mathrm{cir}}}$, which can be rewritten

[TABLE]

which leads to

[TABLE]

It is worth noting that the low-eccentricity approximation allows to cover the range of eccentricities where the criterion (79) is not applicable, since both boundaries depend on the same power of $\varepsilon$.

We plot in Figure 3 the overlap criteria (75) for $\varepsilon=10^{-6}$, the two approximations (77) and (79) from (Deck et al., 2013), as well as the collision condition used in (Laskar & Petit, 2017) approximated for $\alpha\rightarrow 1$,

[TABLE]

We also plot the first MMR islands in order to show the agreement between the proposed criterion and the actual intersections. We see that for high eccentricities, and large $1-\alpha$, the system can verify the MMR overlap stability criterion while allowing for collision between the planets. For small $\alpha$, the MMR overlap criterion alone cannot account for the stability of the system.

4 Critical AMD and MMR

4.1 Critical AMD in a context of resonance overlap

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

The critical AMD was determined thanks to the limit collision condition

[TABLE]

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

We consider a pair as AMD-stable if no distribution of AMD between the two planets allows the overlap of MMR. A first remark is that no pair can be considered as AMD-stable if $\alpha>\alpha_{\mathrm{cir}}$, because in this case, even the circular orbits lead to MMR overlap. Let us write the criterion (75) as a function of $\alpha$ and $\varepsilon$;

[TABLE]

where

[TABLE]

$\sqrt{c_{\mathrm{min}}}$ depends on $\Delta\varpi$ and has a maximum for $\Delta\varpi=\pi$. Since the variation of $\Delta\varpi$ does not affect the AMD of the system, we fix $\Delta\varpi=\pi$ since it is the least-favorable configuration. Therefore we have

[TABLE]

We define the relative AMD of a pair of planets ${\mathscr{C}}$ and express it as a function of the variables $c_{i}$

[TABLE]

The critical AMD $C_{c}^{(1)}$ associated to the overlap criterion (75) can be defined as the smallest value of relative AMD such that the conditions

[TABLE]

are verified by any couple $(c_{1},c_{2})$. We represent this configuration in Figure 4.

As in (Laskar & Petit, 2017), the critical AMD is obtained through Lagrange multipliers

[TABLE]

The tangency condition gives a relation between $c_{1}$ and $c_{2}$,

[TABLE]

Replacing $c_{2}$ in relation (92) gives the critical expression of $c_{1}$ and we immediately obtain the expression of $c_{2}$

[TABLE]

The value of $C_{c}^{(1)}$ is obtained by injecting the critical values $c_{c,1}$ and $c_{c,2}$ into the expression of ${\mathscr{C}}$

[TABLE]

4.2 Comparison with the collision criterion

It is then natural to compare the critical AMD $C_{c}^{(1)}$ to the critical AMD $C_{c}$ (denoted hereafter by $C_{c}^{(0)}$) derived from the collision condition (87). If $\alpha>\alpha_{\mathrm{cir}}$, the circular overlap criterion implies that $C_{c}^{(1)}=0$ and therefore $C_{c}^{(1)}$ should be preferred to the previous criterion $C_{c}^{(0)}$. However, $C_{c}^{(1)}$ was obtained thanks to the assumption that $\alpha$ was close to 1. Particularly, it makes no sense to talk about first-order MMR overlap for $\alpha<0.63$ which corresponds to the center of the MMR 2:1. Therefore, the collision criterion should be used for small $\alpha$. We need then to find $\alpha_{R}$ such that for $\alpha<\alpha_{R}$, we should use the critical AMD $C_{c}^{(0)}$. Since we are close to 1, we use a development of $C_{c}^{(0)}$ presented in (Laskar & Petit, 2017), and similarly, only keep the leading terms in $1-\alpha$ in $C_{c}^{(1)}$. The two expressions are

[TABLE]

We observe that for $\alpha$ close to 1, the two expression have the same dependence on $\gamma$, therefore, $\alpha_{R}$ depends solely on $\varepsilon.$ Simplifying $C_{c}^{(0)}=C_{c}^{(1)}$ gives $\alpha_{R}$ as a solution of the polynomial equation in $(1-\alpha)$;

[TABLE]

While an exact analytical solution cannot be provided, a development in powers of $\varepsilon$ gives the following expression

[TABLE]

It should be remarked that the first term can be directly obtained using Deck’s high-eccentricity approximation.

In Figure 5 we plot $\alpha_{R}$ and $\alpha_{\mathrm{cir}}$ and indicate which criterion is used in the areas delimited by the curves. We specifically represented the region $\alpha>\alpha_{\mathrm{cir}}$ because we cannot treat this region in a similar manner to the remaining region since comparing the relative AMD ${\mathscr{C}}$ to $C_{c}^{(1)}$ does not provide any information. We see that the curve $\alpha_{R}$ is not exactly at the limit where $C_{c}^{(0)}=C_{c}^{(1)}$ for higher $\varepsilon$ due to the development of the critical AMDs for $\alpha\rightarrow 1$. We study the influence of $\gamma$ on the difference between $\alpha_{R}$ and the actual limit in Appendix D

For stability analysis, we need to choose the smallest of the two critical AMD. For $\alpha<\alpha_{R}$, the collisional criterion is better and the MMR overlap criterion is used for $\alpha>\alpha_{R}$. We thus define a piece-wise global critical AMD represented in Figure 6

[TABLE]

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

In (Laskar & Petit, 2017), we proposed a classification of the planetary systems based on their AMD-stability. A system is considered as AMD-stable if every adjacent pair of planets is AMD-stable. A pair is considered as AMD-stable if its AMD-stability coefficient

[TABLE]

where $C$ is the total AMD of the system, $\Lambda^{\prime}$ is the circular momentum of the outer planet and $C_{c}^{(0)}$ is the critical AMD derived from the collision condition. A similar AMD-coefficient can be defined using the global critical AMD defined in (101) instead of the collisional critical AMD $C_{c}^{(0)}$. Let us note $\beta^{(\mathrm{MMR})}$, the AMD-stability coefficient associated to the critical AMD (101).

We can first observe that $\beta^{(\mathrm{MMR})}$ is not defined for $\alpha>\alpha_{\mathrm{cir}}$. Indeed, the conservation of the AMD cannot be guaranteed for orbits experiencing short-term chaos.

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

5.1 Sample and methodology

We first briefly recall the methodology used in (Laskar & Petit, 2017); to which we refer the reader for full details. We compute the AMD-stability coefficients for the systems taken from the Extrasolar Planets Encyclopædia222http://exoplanet.eu/ with known periods, planet masses, eccentricities, and stellar mass. For each pair of adjacent planets, $\varepsilon$ was computed using the expression

[TABLE]

where $m_{1}$ and $m_{2}$ are the two planet masses and $m_{0}$, the star mass. The semi-major axis ratio was derived from the period ratio and Kepler third law in order to reduce the uncertainty.

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

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

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

In (Laskar & Petit, 2017), we studied 131 systems but we did not find the stellar mass for 4 of these systems. They were, as a consequence, excluded from this study. Moreover, the computation of $\varepsilon$ for the pairs of planets of the 127 remaining systems of the sample led in some cases to high planet-to-star mass ratios. We decide to exclude the systems such that $\alpha_{\mathrm{cir}}$ was smaller than the center of the resonance 2:1. We thus discard systems such that a pair of planets has

[TABLE]

As a result, we only consider in this study 111 systems that meet the above requirements.

A pair is considered stable if the 1-$\sigma$ confidence interval (84% of the simulated systems) of the AMD-stability coefficient $\beta^{(\mathrm{MMR})}$ is below 1. A system is stable if all adjacent pairs are stable.

5.2 Results

Figure 7 shows the planet pairs of the considered systems in a plane $\alpha$-$\varepsilon$. The color associated to each point is the AMD-stability coefficient of the pair. The values chosen for the plot correspond for all quantities to the median. We remark that very few systems are concerned by the change of the critical AMD, indeed, only eight systems333 It should be noted that for one of the systems, the MMR overlap criterion was preferred in 16% of the Monte Carlo simulations. have a pair of planets such that $C_{c}^{(1)}<C_{c}^{(0)}$. The 111 considered systems contain 162 planet pairs plotted in Figure 7. This means that less than 5% of the pairs are in a configuration leading to MMR overlap.

We plot in Figure 8, the architecture of these eight systems and give in Table E.1 the values of the AMD-stability coefficients. For each of these systems, the pair verifying the MMR overlap criterion was already considered AMD-unstable by the criterion based on the collision.

In order to show this, we plot in Figure 9 the AMD-stability coefficients computed with both critical AMD. We see that the pairs affected by the change of criterion were already considered AMD-unstable in the purely secular dynamics. However, while those pairs have a collisional AMD-coefficient $\beta$ between 1 and 10, the global AMD-stability coefficient is increased by roughly an order of magnitude for the four pairs with $\alpha$ between $\alpha_{R}$ and $\alpha_{\mathrm{cir}}$. The AMD-coefficient is not defined for the three pairs verifying the circular MMR overlap criterion. The pair HD 47366 b/c does not see a significant change of its AMD-stability coefficient due to the small number of cases where $C_{c}^{(0)}>C_{c}^{(1)}$.

We identify three systems, HD 200964, HD204313 and HD 5319, that satisfy the circular overlapping criterion. As already explained in (Laskar & Petit, 2017), AMD-unstable planetary systems may not be dynamically unstable. However, it should be noted that the period ratios of the AMD-unstable planet pairs are very close to particular MMR.

Indeed, we have

[TABLE]

The AMD-instability of those systems strongly suggests that they are indeed into a resonance which stabilizes their dynamics.

6 Conclusions

As shown in Laskar & Petit (2017), the notion of AMD-stability is a powerful tool to characterize the stability of planetary systems. In this framework, the dynamics of a system is reduced to the AMD transfers allowed by the secular evolution.

However, we need to ensure that the system dynamics can be averaged over its mean motions. While a system can remain stable and the AMD or semi-major axis can be averaged over timescales longer than the libration period in presence of MMR, the system stability and particularly the conservation of the AMD is no longer guaranteed if the system experiences MMR overlap. In this paper, we use the MMR overlap criterion as a condition to delimit the zone of the phase space where the dynamics can be considered as secular.

We refine the criteria proposed by (Wisdom, 1980; Mustill & Wyatt, 2012; Deck et al., 2013) and demonstrate that it is possible to obtain a global expression (75), valid for all cases. The previous circular (77) and eccentric (79) criteria an then be derived from (75) as particular approximations. Moreover, we show that expression (75) can be used to directly take into account the first-order MMR in the notion of AMD-stability.

With this work on first-order MMR, we improve the AMD-stability definition by addressing the problem of the minimal distance between close orbits. For semi-major axis ratios $\alpha$ above a given threshold $\alpha_{\mathrm{cir}}$ (77), that is, $\alpha_{\mathrm{cir}}<\alpha<1$, the system is considered unstable whichever value the AMD may take given that even two circular orbits satisfy the MMR overlap criterion. At wider separations, circular orbits are stable but as eccentricities increase two outcomes may happen: Either the system enters a region of MMR overlap or the collision condition is reached. The system is said to be AMD-unstable as soon as any of these conditions is reached. Above a second threshold, $\alpha_{R}<\alpha<\alpha_{\mathrm{cir}}$ (Eq. 100) the AMD-stability is governed by MMR overlap while for wider separations ($\alpha<\alpha_{R}$) we retrieve the critical AMD defined in (Laskar & Petit, 2017) which only depends on the collision condition.

In order to improve the AMD-stability definition for the collision region, we could even take into account the non-secular dynamics induced by higher-order MMR and close-encounter consequences on the AMD. To study this requires more elaborated analytical considerations than those presented here that are restricted to the first-order MMR; this will be the goal of future work.

We show in Section 5 that very few systems satisfy the circular MMR overlap criterion. Moreover, the presence of systems satisfying this criterion strongly suggests that they are protected by a particular MMR. In this case, the AMD-instability is a simple tool suggesting unobvious dynamical properties.

Appendix A Expression of the first-order resonant Hamiltonian

We use the method proposed in (Laskar 1991) and (Laskar & Robutel 1995) to determine the expression of the planetary perturbation $\hat{\mathcal{H}}_{1}$. $\hat{\mathcal{H}}_{1}$ can be decomposed into a part from the gravitational potential between planets $\hat{U}_{1}$ and a kinetic part $\hat{T}_{1}$ as

[TABLE]

with

[TABLE]

The difficulty comes from the development of $a_{2}/\Delta_{12}$ and its expression in terms of Poincaré variables. We note $S$, the angle between ${\bf u}_{1}$ and ${\bf u}_{2}$. We have

[TABLE]

Let us denote $\rho=u_{1}/u_{2}$, $a_{2}/\Delta_{12}$ can be rewritten

[TABLE]

where we denote

[TABLE]

$V$ is at least of order one in eccentricity. We can therefore develop (112) for small $V$. We only keep the terms of first order in eccentricity,

[TABLE]

The well-known development of the circular coplanar motion $A$ gives (e.g., Poincaré 1905)

[TABLE]

where $b^{(k)}_{s}(\alpha)$ are the Laplace coefficients (25).

Because of the averaging over the non-resonant fast angles, the non-vanishing terms have a dependence on $\lambda_{i}$ of the form $j\left((p+1)\lambda_{2}-p\lambda_{1}\right)$. Since we only keep the terms of first order in eccentricity, the d’Alembert’s rule (7) imposes $j=\pm 1$. Let us compute the first-order development of $a_{2}/u_{2}$ and $V$ in terms of Poincaré variables and combine these expressions with $A^{-1/2}$ and $A^{-3/2}$ in order to select the non-vanishing terms.

Let us denote $z_{i}={\rm e}^{\mathbf{i}\lambda_{i}}$ and $z=z_{1}\bar{z}_{2}={\rm e}^{\mathbf{i}(\lambda_{1}-\lambda_{2})}$. The researched terms are of the form

[TABLE]

Let us denote

[TABLE]

the first term in the development (117) gives

[TABLE]

The contributing term has for expression

[TABLE]

where $\mathbf{\mathfrak{X}}_{i}=\hat{\mathfrak{\scriptstyle X}}_{i}\sqrt{2/\hat{\Lambda}_{i}}=X_{i}\mathbf{e}^{\mathbf{i}((p+1)\lambda_{2}-p\lambda_{1})}$.

For the computation of the second term of (117), the only contribution comes from $V$ since $a_{2}/u_{2}\sim 1$. We define

[TABLE]

$V$ can be expressed as a function of $z,\bar{z},U$ and $\bar{U}$. Indeed we have

[TABLE]

and

[TABLE]

where $\mathrm{O}(e^{2})$ corresponds to terms of total degree in eccentricities of at least 2. We deduce from these two last expressions that

[TABLE]

We can therefore write444 In Laskar & Robutel (1995) the first-order expression of $V$ is written $W_{1}=(UZ+\bar{U}\bar{Z})$ instead of ${W_{1}=(UZ+\bar{U}\bar{Z})/2}$. This misprint in equation (47) of (Laskar & Robutel 1995) is transmitted as well in equation (51). It has no consequences in the results of the paper.

[TABLE]

where $Z=\alpha(3\bar{z}-2\alpha-z)$. With this expression of $V$, it is easy to gather the corresponding terms and the second term in the development (117) gives the contributing term

[TABLE]

After gathering the terms (123,130), we can give the expression of the resonant Hamiltonian

[TABLE]

where

[TABLE]

with $\gamma=m_{1}/m_{2}$, and

[TABLE]

The kinetic part $\hat{T}_{1}$ has no contribution to the averaged resonant Hamiltonian for $p>1$. Indeed, as explained above, due to the d’Alembert rule, the first-order terms must have an angular dependence of the form $j(-p\lambda_{1}+(p+1)\lambda_{2})$. At the first order in $\varepsilon$, such a term can only be present in the development of the inner product ${\bf\tilde{u}}_{1}\cdot{\bf\tilde{u}}_{2}$. At the first order in eccentricities, we have (Laskar & Robutel 1995)

[TABLE]

where $\omega_{j}$ is the true longitude of the planet $j$. The only term with the good angular dependence comes from $\Re\mathbf{e}^{\mathbf{i}(\omega_{1}-\omega_{2})}$ since the other first-order terms only depend on one mean longitude. The development of $\mathbf{e}^{\mathbf{i}(\omega_{1}-\omega_{2})}$ at the first order in eccentricities gives

[TABLE]

Thus for $p>1$, $\hat{T}_{1}$ has no contribution to the averaged Hamiltonian, and for $p=1$ we have

[TABLE]

A.1 Asymptotic expression of the resonant coefficients

We present the method we used to obtain the analytic development of the coefficients $r_{1}$ and $r_{2}$ defined in equations (135) and (136). Using the expression of $b^{(k)}_{s}(\alpha)$, we have

[TABLE]

We can rewrite this expression

[TABLE]

We make the change of variable $\phi=(1-\alpha)u$ in the integrals. Factoring $(1-\alpha)^{3}$, the denominators in the integrals can be developed for $\alpha\rightarrow 1$

[TABLE]

Using the relation $\alpha_{0}=(p/(p+1))^{2/3}$, the numerators can be developed

[TABLE]

Therefore, we deduce the equivalent of $r_{1}$ for $p\rightarrow+\infty$

[TABLE]

where $K_{\nu}(x)$ is the modified Bessel function of the second kind. Similarly, we have $r_{2}\sim-r_{1}$ since the additional term is of lower order in $p$.

We can obtain the constant term of the development by using the second order expression of $\alpha_{0}$ and developing the integrand to the next order in $(1-\alpha)$. We give here the numerical expressions of the two developments

[TABLE]

Appendix B Development of the Keplerian part

We show here that the first order in $(C-\Delta G)$ of the Keplerian part vanishes and give the details of the computation for the second order. The Keplerian part can be written

[TABLE]

Therefore, the first order in $C-\Delta G$ has for expression

[TABLE]

since we have

[TABLE]

The second-order term has for coefficient

[TABLE]

Appendix C Width of the resonance island

We detail in this Appendix the computation of the resonance island’s width (see also Ferraz-Mello 2007, Appendix C).

C.1 Coefficients-roots relations

We first explain how the width of the resonance can be related to the position of the saddle point on the $X$-axis. The resonant island has a maximal width on the $X$-axis. Therefore we need to compute the expression of the intersections of the separatrices with the $X$-axis.

Let us note $\mathcal{H}_{3}$, the energy at the saddle point $(X_{3},0)$. Since the energy of the separatrices is $\mathcal{H}_{3}$ as well, the two intersections of the separatrices with the $X$-axis are the solution of the equation

[TABLE]

This equation has three solutions $X_{1}^{*},X_{2}^{*}$, and $X_{3}$ which has a multiplicity of 2. We can therefore rewrite the equation as

[TABLE]

We detail here the relations between the coefficients and the roots of the polynomial equation (154). We have

[TABLE]

From relation 155, we have directly $X_{1}^{*}+X_{2}^{*}=-2X_{3}$, and since

[TABLE]

we can express $(X_{1}^{*}-X_{2}^{*})^{2}$ as a function of $X_{3}$ thanks to the relations (156) and (157)

[TABLE]

We thus deduce the expressions of $X_{1}^{*}$ and $X_{2}^{*}$ as functions of $X_{3}$

[TABLE]

As explained in section 3.1, we obtain the width of the resonance in terms of variation of $\alpha$ as a function of $X_{3}$ (equation (61)). We can use this expression to obtain the width of the resonance for particular cases detailed in the following subsections.

C.2 Width for initially circular orbits

In the case of initially circular orbits, the minimal AMD to enter the resonance is 0. For ${C_{\mathrm{min}}}=0$, the equation (63) gives $X_{3}=2^{2/3}$ as a solution and we have

[TABLE]

We find here the same width of resonance as (Deck et al. 2013).

C.3 Width for highly eccentric orbits

If we consider a system with ${C_{\mathrm{min}}}\gg\chi^{2/3}$, our formalism gives us the result first proposed by Mustill & Wyatt (2012) and improved by Deck et al. (2013) for eccentric orbits. In this case, we can inject the approximation (67) of $X_{3}$ in the expression (61) of $\delta\alpha$ and obtain

[TABLE]

This result is also similar to Deck’s one, using $\sqrt{c_{\mathrm{min}}}$ instead of $\sigma$ (Deck et al. 2013, equation (25)).

C.4 Width for low eccentric orbits

For ${C_{\mathrm{min}}}\ll\chi^{2/3}$, we propose here a new expression of the width of resonance thanks to the expression (68). This expression is an extension of the circular result presented above (162). Let us develop $\sqrt{X_{3}}$ for ${C_{\mathrm{min}}}\ll\chi^{2/3}$

[TABLE]

Therefore for low-eccentricity systems, we have

[TABLE]

where $\delta\alpha_{c}$ is the width of the resonance for initially circular orbits defined in (162).

Appendix D Influence of $\gamma$ on the limit $\alpha_{R}$

As can be seen in Figure 5, the solution $\alpha_{R}$ of equation (99) is not the exact limit where the collision and the MMR criteria are equal. Indeed, equation (99) is obtained after the development of $C_{c}^{(0)}$ and $C_{c}^{(1)}$ for $\alpha$ close to 1. Since at first order, both expressions have the same dependence on $\gamma$, $\alpha_{R}$ does not depend on $\gamma$. In order to study the dependence on $\gamma$ of the limit $\alpha_{\mathrm{lim}}$ where $C_{c}^{(0)}=C_{c}^{(1)}$, we plot in Figure (10), for different values of $\varepsilon$, the quantity

[TABLE]

which gives the error made when approximating $\alpha_{\mathrm{lim}}$ by $\alpha_{R}$. We see that all the curves have the same shape with an amplitude increasing with $\varepsilon$. For high $\gamma$, $\alpha_{R}$ is very accurate even for the greatest values of $\varepsilon$. Moreover, the error is maximum for very small $\gamma$ and always within a few percent.

The amplitude of the error scales with ${1-\alpha_{R}\propto\varepsilon^{1/4}}$ as we can see in the Figure 11. We plot in this Figure 11 the quantity $\delta\alpha_{R}/\varepsilon^{1/4}$; we see that the curves are almost similar, particularly for the smaller values of $\varepsilon$.

Appendix E AMD-stability coefficients of the system affected by the MMR overlap criterion

We report in Table E.1 the AMD-stability coefficients of the systems where more than 5% of the Monte Carlo realizations were affected by the change of critical AMD. Apart for the system HD 47366 where 16% of the simulations used the new criterion, the seven other systems used the critical AMD $C_{c}^{(1)}$ for almost all the realizations. For HD 204313, only the pair (b/d) is affected.

In Table E.1, $\sqrt{\langle e^{2}\rangle}$ corresponds to the mean value of the squared eccentricity computed as explained in section (5.1).