© ESO 2018

1 Introduction

The chaotic nature of planetary dynamics has made the development of stability criteria necessary because the complete study of individual systems for the lifetime of the central star would require too much computational power. Moreover, in the context of exoplanet dynamics, the orbital parameters are often not known with great precision, making it impossible to conduct a precise dynamical study. For multiplanetary systems, the best solutions yet are empirical stability criteria based on minimum spacing between planets obtained from numerical simulations (Chambers et al. 1996; Pu & Wu 2015). For tightly packed systems, Tamayo et al. (2016) have suggested a solution based on short integrations and a machine-learning algorithm.

Another approach consists in the study of the angular momentum deficit (AMD; see Laskar 1997, 2000). The AMD is a weighted sum of the eccentricities and mutual inclinations of planets and can be interpreted as a dynamical temperature of the planetary system. The AMD is conserved at all orders of averaging over the mean motions. In Laskar & Petit (2017), it was shown that if the AMD is small enough, collisions are impossible. We can define a sufficient stability condition from the secular conservation of the AMD, i.e., the AMD-stability. Systems that are AMD-stable are long-lived, whereas a more in-depth dynamical study is necessary for AMD-unstable systems.

The initial AMD-stability definition is based on the secular approximation. In Petit et al. (2017), the criterion was slightly modified to exclude systems experiencing short-term chaos due to first-order mean motion resonances (MMR) overlap. The MMR overlap is the main source of chaos and instability in planetary dynamics. Based on Chirikov (1979) resonance overlap criterion, Wisdom (1980) has derived an analytical stability criterion for the two-planet systems that has since been widely used. Wisdom’s criterion is obtained for circular and coplanar planets. It was improved to take into account moderate eccentricities (Mustill & Wyatt 2012; Deck et al. 2013). We have shown in Petit et al. (2017) that this MMR overlap criterion can be expressed solely as a function of the AMD, semimajor axes, and masses for almost coplanar and low eccentricity systems.

The MMR overlap criterion gives a clear limit between regular and chaotic orbits for small eccentricities. However, this criterion is based on first-order expansions in the planets-to-star mass ratio, the spacing between planets, eccentricities, and inclinations. As a result, for wider orbit separations, the secular collision criterion from Laskar & Petit (2017) remains a better limit (Petit et al. 2017). Moreover, the MMR overlap criterion only takes into account the interaction between a couple of planets, making it really accurate only for two-planet systems.

In the case of two-planet systems, the topology of the phase space gives a far simpler criterion of stability, that is the Hill stability. Based on a work by Sundman (1912) on the moment of inertia in the three-body problem, Marchal & Saari (1975) noted the existence of forbidden zones in the configuration space. In Marchal & Bozis (1982), they extended the notion of Hill stability to the general three-body problem and showed that some systems can forbid close encounters between the outer body and any of the inner bodies. In particular, the Hill stability ensures that collisions between the outer body and the close binary (in a two-planet system, between the outer planet and the inner planet or the star) are impossible for infinite time.

The results of Marchal & Saari (1975) and Marchal & Bozis (1982) have many applications outside of the Hill stability. It is possible to cite a sufficient condition for the ejection of a body from the system (Marchal et al. 1984a,b) or the determination of the limit of a triple close approach for bounded orbits (Laskar & Marchal 1984).

Marchal & Bozis (1982) presented the planetary problem (one body with a much larger mass) as a particular case, but the result was mainly popularized by Gladman (1993) who introduced a minimal spacing for initially circular and coplanar systems. Gladman also proposed some criteria for eccentric orbits in some particular configurations of masses. The result by Gladman was refined to cover other situations as, for example, the case of inclined orbits (Veras & Armitage 2004). Georgakarakos (2008) provided a review of stability criteria for hierarchical three-body problems.

The Hill stability is consistent with numerical integrations, where a sharp transition between Hill stable and Hill unstable systems is often observed (Gladman 1993; Barnes & Greenberg 2006; Deck et al. 2013). Moreover, Deck et al. (2013) analyzed the differences between the MMR overlap and Hill stability criteria. They remarked that there exists an area where orbits are chaotic because of the overlap of MMR but long-lived owing to Hill stability.

In this paper, we show that Marchal & Bozis (1982) Hill stability criterion fits extremely well in the AMD-stability framework. In Sect. 2, we derive a criterion for Hill stability solely expressed as a function of the total AMD, semimajor axes, and masses. Our criterion does not need any expansion in the spacing, eccentricities, or inclinations of the orbits and admits all previous Hill criteria as particular approximations. In order to do so, we follow the reasoning proposed by Marchal & Bozis (1982).

We then compare the Hill stability criterion with the AMD-stability criteria proposed in Laskar & Petit (2017) and Petit et al. (2017). In the last section, we carry numerical integrations of two-planet systems over a large part of the phase space. We show that the only parameters of importance are the initial AMD and semimajor axis ratio.

2 Hill stability in the three-body problem

2.1 Generalized Hill curves

Let us use the formalism proposed in Marchal & Saari (1975) and Marchal & Bozis (1982) for the definition of the Hill regions for the general three-body problem. We mainly consider the planetary case and therefore slightly adapt their notations to this particular problem. Let us consider two planets of masses m₁, m₂ orbiting a star of mass m₀. Let $\mathcal{G}$ be the constant of gravitation, $\mu =\mathcal{G}{m}_0$, $$\epsilon =\frac{m_1+m_2}{m_0},$$(1)

the planet mass to star mass ratio and $$\gamma =\frac{m_1}{m_2}.$$(2)

As in Laskar & Petit (2017), we use the heliocentric canonical coordinates ${({\stackrel{\sim }{r}}_j,{\text{r}}_j)}_{j=1,2}$. In those coordinates, the Hamiltonian is $$\mathcal{H}={\sum }_{j=1}^2\left(\frac{1}{2}\frac{{\Vert {\stackrel{}{r\sim }}_j\Vert }^2}{m_j}-\frac{\mu m_j}{r_j}\right)+\frac{1}{2}\frac{{\Vert {\stackrel{}{r\sim }}_1+{\stackrel{}{r\sim }}_2\Vert }^2}{m_0}-\frac{\mathcal{G}{m}_1{m}_2}{r_{12}},$$(3)

where r₁₂ = ∥r₁ −r₂∥. We denote G the total angular momentum of the system: $$G={\sum }_{j=1}^2{\text{r}}_j\times {\stackrel{}{r\sim }}_j,$$(4)

and G its norm. From now, we assume G to be aligned along the z-axis. We also define $${\Lambda }_j=m_j\sqrt{\mu a_j},$$(5)

where a_j is the semimajor axis of the jth planet. We also note e_j, the eccentricity, and i_j, the inclination of the orbital plane with the horizontal plane. The AMD, C (Laskar 1997, 2000) is expressed $$C={\sum }_{j=1}^2{\Lambda }_j\left(1-\sqrt{1-e_j^2}\cos (i_j)\right).$$(6)

Following Marchal & Bozis (1982), we define a generalized semimajor axis: $$a=-\frac{\mathcal{G}{M}^{*}}{2\mathcal{H}},$$(7)

and a generalized semilatus rectum: $$p=\frac{m_0(1+\epsilon )}{\mathcal{G}{M}^{*2}}{G}^2,$$(8)

where $$M^{*}=m_0{m}_1+m_0{m}_2+m_1{m}_2=m_0^2\epsilon \left(1+\epsilon \frac{\gamma }{{(1+\gamma )}^2}\right).$$(9)

The values a and p are the two length units that can be built from the first integrals $\mathcal{H}$ and G.

We finally define two variable lengths. First, ρ the mean quadratic distance $$M^{*}{\rho }^2=m_0{m}_1{r}_1^2+m_0{m}_2{r}_2^2+m_1{m}_2{r}_{12}^2,$$(10)

that is proportional to the moment of inertia I computed in the center of mass frame (Marchal & Bozis 1982): $$I=\frac{1}{2}\frac{M^{*}{\rho }^2}{m_0+m_1+m_2}.$$(11)

We also define ν, the mean harmonic distance: $$\frac{M^{*}}{\nu }=\frac{m_0{m}_1}{r_1}+\frac{m_0{m}_2}{r_2}+\frac{m_1{m}_2}{r_{12}},$$(12)

which is proportional to the potential energy: $$U=-\frac{\mathcal{G}{M}^{*}}{\nu }.$$(13)

For a system with given $\mathcal{H}$ and G, some configurations of the planets are forbidden. Indeed, the value of the ratio ρ∕ν is constrained by the inequality (Marchal & Saari 1975): $$\frac{\rho }{\nu }\ge \frac{\rho }{2a}+\frac{p}{2\rho },$$(14)

derived from Sundman’s inequality (Sundman 1912). Moreover, if the system has a negative energy, the right-hand side of Eq. (14) has a minimum value obtained for $\rho =\sqrt{ap}$. We therefore have the inequality (Marchal & Bozis 1982): $$\frac{{\rho }^2}{{\nu }^2}\ge \frac{p}{a}=-\frac{2m_0(1+\epsilon )}{{\mathcal{G}}^2{M}^{*3}}\mathcal{H}{G}^2.$$(15)

If p∕a is high enough, the inequality (15) makes some regions of the phase space inaccessible. In this case, we can ensure that certain initial conditions forbid collisions between the two planets for all times.

Let us study the values of the function ${(\rho /\nu )}^2$. Since this ratio only depends on the ratios of mutual distances, we can always place ourselves in the planegenerated by the three bodies. We can also choose to place the first planet on the x-axis and normalize the lengths by r₁. Let us call this plane $\mathcal{P}$ and note (x, y) the coordinates of the second planet. In the plane $\mathcal{P}$ (see Fig. 1), the star S is at the origin, the first planet P₁ is situated at the point (1, 0). For the planar restricted three-body problem, this reduction is equivalent to study the dynamics in the corotating frame. We note $$R(x,y)={\left(\frac{\rho }{\nu }\right)}^2.$$(16)

The shape of the function R in $\mathcal{P}$ only depends on the mass distribution, for example, the two ratios ε and γ. In Fig. 1, we can see level curves of the function R plotted in $\mathcal{P}$ for ε = 10⁻³ and γ = 1. The value R is minimal and equal to 1 at the two Lagrange points L₄ and L₅ and goes to + ∞ for ∥r₂∥→ 0, ∥r₂ −r₁∥→ 0 or ∥r₂ ∥ → +∞. The function has three saddle nodes at the Lagrange points L₁, L₂, and L₃. We give the method to compute the position of the Lagrange points in Appendix A. Let us denote x₁, x₂, and x₃ their abscissa, respectively. At the lowest order in ε, we have $$\begin{array}{ll}{x}_1\hfill & =1-{\left(\frac{\epsilon }{3}\right)}^{1/3}+\text{O}({\epsilon }^{2/3}),\hfill \\ x_2\hfill & =1+{\left(\frac{\epsilon }{3}\right)}^{1/3}+\text{O}({\epsilon }^{2/3}),\hfill \\ x_3\hfill & =-1+\frac{7}{12}\frac{\gamma -1}{\gamma +1}\epsilon +\text{O}({\epsilon }^2).\hfill \end{array}$$(17)

At the lowest order in ε, the value of R at the Lagrange points are as follows: $$\begin{array}{ll}R(L_1)\hfill & =1+3^{4/3}{\epsilon }^{2/3}\gamma /{(\gamma +1)}^2+\text{O}(\epsilon ),\hfill \\ R(L_2)\hfill & =1+3^{4/3}{\epsilon }^{2/3}\gamma /{(\gamma +1)}^2+\text{O}(\epsilon ),\hfill \\ R(L_3)\hfill & =1+2\epsilon \gamma /{(\gamma +1)}^2+\text{O}({\epsilon }^2).\hfill \end{array}$$(18)

The values R(L₁) and R(L₂) have the same first-order term but differ in the expansion of higher order. Indeed, if m₀ ≥ m₁ ≥ m₂, we have R(L₁) ≥ R(L₂) (Marchal & Bozis 1982). From now on, let us assume R(L₁) ≥ R(L₂) (if not, we can just substitute R(L₂) to R(L₁) in further equations).

For p∕a ≥ R(L₁), the accessible domain is split into three parts: the Hill sphere of the star S, which is around the origin; the first planet Hill sphere¹ $S_{{\text{H}}_{\text{1}}}$ (in green in Fig. 1); and the outer region. In this case, if the second planet is not initially inside $S_{{\text{H}}_{\text{1}}}$, it will never be able to enter this region. However, if P₂ is in the outer region, the Hill stability cannot constrain the possibility of ejection. Similarly, if P₂ is closer to the star (in the inner region), a collision with the star is still possible.

The study of the function R and the inequality (15) gives a noncollision criterion for an infinite time. Marchal and Bozis called it the Hill stability.

$$\frac{p}{a}>R(L_1),$$(19)

where a is defined in (7), p in (8), and R in (16).

From this inequality, Gladman (1993) obtained criteria for initially circular orbits and for two particular cases of eccentric orbits: the case of equal masses and small eccentricities and the case of equal masses and large, but equal eccentricities.

While Gladman’s Hill stability criterion for initially circular orbits is useful, the eccentric criteria are too particular to be used in the context of a generic system. It is, however, possible to obtain a very general Hill stability criterion using the AMD to take into account the eccentricities and inclinations of the orbits.

Fig. 1 Some levels of the function R defined in (16) for ε = 10−3 and γ = 1 in the $\mathcal{P}$ plane. The two red points correspond to the Lagrange points L1 and L2 and the three orange points to L3, L4, and L5. The orange-filled area corresponds to the region where R(x, y) < R(L1). The star S is at the origin, the planet P1 at (1,0), and (x,y) are the coordinates of the second planet. The green region $S_{{\text{H}}_{\text{1}}}$represents the Hill sphere of the first planet.

2.2 AMD condition for Hill stability

The total energy of the system can be written $$\mathcal{H}=-\frac{m_2^3{\mu }^2}{2{\Lambda }_2^2}\left(\frac{\gamma }{\alpha }+1+h_1\right),$$(20)

where α = a₁∕a₂ and $$h_1=-\frac{2{\Lambda }_2^2}{m_2^3{\mu }^2}\left(\frac{1}{2}\frac{{\Vert {\stackrel{}{r\sim }}_1+{\stackrel{}{r\sim }}_2\Vert }^2}{m_0}-\frac{\mathcal{G}{m}_1{m}_2}{r_{12}}\right).$$(21)

From now on, we assume that initially α ≤ 1 (if not we can just renumber the two planets). Similarly, the angular momentum can be rewritten: $$G={\Lambda }_2\left(\gamma \sqrt{\alpha }+1-\mathcal{C}\right),$$(22)

where $$\mathcal{C}=\frac{\mathcal{C}}{{\Lambda }_2}=\gamma \sqrt{\alpha }\left(1-\sqrt{1-e_1^2}\cos i_1\right)+1-\sqrt{1-e_2^2}\cos i_2$$(23)

is the relative AMD (Laskar & Petit 2017). Combining Eqs. (15), (20), and (22), we obtain $$\frac{p}{a}=\frac{(1+\epsilon )}{{(\gamma +1)}^3{(1+\epsilon \gamma /{(\gamma +1)}^2)}^3}\left(\frac{\gamma }{\alpha }+1+h_1\right){\left(\gamma \sqrt{\alpha }+1-\mathcal{C}\right)}^2.$$(24)

The Hill stability criterion (19) can be rewritten without any approximation as a condition on $\mathcal{C}$ and we have the following formulation of the Hill stability.

$$\mathcal{C}<C_c^{\text{Ex}}=\gamma \sqrt{\alpha }+1-{(\gamma +1)}^{3/2}\sqrt{\frac{R(L_1){(1+\epsilon \gamma /{(\gamma +1)}^2)}^3}{(1+\epsilon )\left(\gamma /\alpha +1+h_1\right)}},$$(25)

where $\mathcal{C}$ is the relative AMD (23) and h₁ is the normalized perturbation part (21).

The inequality (25) is equivalent to Proposition 1, but we isolated the contribution from the AMD on the left-hand side. Up to the perturbation term h₁, the right-hand side of Eq. (25) only depends on the masses and the semimajor axis ratio α. If we only keep the terms of leading order in ε in the square root of the right-hand side of Eq. (25), we obtain an expression that depends only on α, ε, and γ.

$$\mathcal{C}<\gamma \sqrt{\alpha }+1-{(1+\gamma )}^{3/2}\sqrt{\frac{\alpha }{\gamma +\alpha }\left(1+\frac{3^{4/3}{\epsilon }^{2/3}\gamma }{{(1+\gamma )}^2}\right)}+\text{O}(\epsilon ).$$(26)

As explained in Appendix B, h₁ is of smaller order in ε and can be neglected if the criterion is verified. We want to stress that the expression (26) is obtained with only an expansion in ε and only depends on α, $\mathcal{C}$, and the masses of the bodies. The term of order ε^2∕3 in Eq. (26) also depends on $\gamma /{(\gamma +1)}^2$, but we show in Appendix B that (26) is still valid for γ ≪ 1 or γ ≫ 1.

2.3 Close planets approximation

Assuming 1 − α ≪ 1, and a smallAMD value, further approximations can be made. At leading order in $\mathcal{C},1-\alpha$, and ε, the inequality (26) becomes $$\mathcal{C}<\frac{3\gamma }{8(\gamma +1)}{(1-\alpha )}^2-\frac{3^{4/3}\gamma }{2(\gamma +1)}{\epsilon }^{2/3}.$$(27)

We can isolate 1 − α in this expression to obtain an approximate minimum spacing for Hill stable systems.

$$\frac{a_2-a_1}{a_2}=1-\alpha >\sqrt{4\times 3^{1/3}{\epsilon }^{2/3}+\frac{8}{3}\frac{\gamma +1}{\gamma }\mathcal{C}}.$$(28)

Gladman’s eccentric criteria can be recovered from Eq. (28) if the AMD is developed under the assumptions (same planet masses, small or large, and equal eccentricities) made in Gladman (1993). However, Eq. (28) is more general as it takes into account mutual inclinations or uneven mass distribution.

In the case of circular orbits, we also get the well-known formula (Gladman 1993): $$1-\alpha >2\times 3^{1/6}{\epsilon }^{1/3}=2.40{\epsilon }^{1/3}.$$(29)

2.4 Comparison of the Hill criteria

We can compare the right-hand side of Eqs. (25)–(27), and Gladman’s circular approximation (29) to test how relevant are the approximations made here. In Fig. 2, we plot the exact expression $C_c^{\text{Ex}}$ (Eq. (25); in green), the expansion in ε (Eq. (26); in orange), the approximation for close planets (Eq. (27); in red), and the minimum spacing for circular orbits (Eq. (29); in blue). We see that the expansion in ε (Eq. (26)) cannot be distinguished from the exact curve (Eq. (25)). In order to better quantify this, we plot in Fig. 3 the maximum difference between the two curves as a function of ε for various values of the mass ratio γ.

We seein Fig. 3 that for the range of ε used in planetary dynamics (typically from 10⁻⁶–10⁻²), the expression (26) developed in ε is accurate even for very uneven planet mass distribution. From now, we use Eq. (26) to define the Hill stability.

Fig. 2 Comparison of the right-hand sides of the inequalities – Eq. (25) (orange), Eq. (26) (green), and Eq. (27) (red) – as a function of α for γ = 1 and ε = 10−5 (upper three curves) or ε = 10−3 (lower three curves). The black curve is the critical AMD for the collision condition (Laskar & Petit 2017). The green and orange curves are on top of each other (see Fig. 3). The critical AMD from MMR overlap$C_{c^{\text{MMR}}}$ and circular criterion of Gladman (1993) are also plotted for comparison.

Fig. 3 Maximum difference between $C_c^{\text{Ex}}$ (Eq. (25)) and $C_c^{\text{H}}$ (Eq. (30)) on the right-hand sides of Eqs. (25) and (26), respectively, as a function of ε for various values of γ = m1∕m2. In Eq. (25), the term h1 is evaluated using the approximation of the kinetic term given in Appendix B.

3 Comparison with the AMD-stability

In Laskar &Petit (2017), the AMD-stability of a system is defined by comparing its AMD to a critical value C_c above which the stability of the system cannot be guaranteed. The AMD-stability is a sufficient condition for long-lived stability of a system. The result from Laskar & Petit (2017) was obtained by finding the minimal AMD needed for a system to allow for collisions in the secular approximation, i.e., assuming the semimajor axes are fixed. We completed, in Petit et al. (2017), the AMD-stability definition by excluding the configurations such that the overlap of first-order MMR occurs. The modified critical AMD is hereafter noted $C_{c^{\text{MMR}}}$.

We can use Eq. (26) to define a critical AMD $C_c^{\text{H}}$ for the Hill stability: $$C_c^{\text{H}}=\gamma \sqrt{\alpha }+1-{(1+\gamma )}^{3/2}\sqrt{\frac{\alpha }{\gamma +\alpha }\left(1+\frac{3^{4/3}{\epsilon }^{2/3}\gamma }{{(1+\gamma )}^2}\right)}.$$(30)

A system isHill stable if its initial relative AMD is smaller than the initial critical AMD $C_c^{\text{H}}$.

We see that for two planets, the Hill stability definition fits extremely well in the AMD-stability framework. We can also compare $C_c^{\text{H}}$ to the previously proposed critical AMD. The collision critical AMD C_c is plotted in Fig. 2 with two values of $C_c^{\text{H}}$ for ε = 10⁻⁵ (resp. 10⁻³). We can see that the Hill stability criterion is stricter ($C_c^{\text{H}}<C_c$) than the collision condition for secular dynamics. It can be easily understood since the Hill stability forbids the planets to approach each other.

Indeed, let us consider a Hill stable system, i.e., such that $\mathcal{C}<C_c^{\text{H}}$. As a result, the two planets cannot approach each other by less than their mutual Hill radius for any variation of semimajor axes and thus also in the secular system. In particular, a configuration such that the two orbits intersect is impossible. Therefore, the system is AMD-stable and $\mathcal{C}<C_c$. Since we are not making any additional hypothesis regarding $\mathcal{C}$, we have $C_c^{\text{H}}<C_c$. The strict inequality comes from the positive minimal distance between the two planets.

As a comparison, we also plot in Fig. 2 the MMR critical AMD $C_{c^{\text{MMR}}}$. For small relative AMD, $\mathcal{C}$, $C_{c^{\text{MMR}}}$, and $C_c^{\text{H}}$ are almost identical.

4 Numerical simulations

The Hill criterion proposed by Marchal & Bozis (1982) has already been tested numerically in particular cases (Gladman 1993; Veras & Armitage 2004; Barnes & Greenberg 2006). In their comparison between the Hill and the overlap of MMR criteria, Deck et al. (2013) noted a sharp transition in the proportion of chaotic orbits at the Hill limit (p∕a) = (p∕a)|_c. It also appears that for small ε and α close to 1, the overlap of MMR criterion provides a better limit for the chaotic region.

We want to test whether the Hill criterion gives a good limit to the chaotic region for wider separations. Moreover, we want our initial conditions to sample homogeneously the phase space. Indeed, the Hill stability criterion studied in this paper only depends on few quantities, the relative AMD $\mathcal{C}$, and the ratio of semimajor axis and not the angles or the actual distribution of the AMD between the degrees of freedom of eccentricities or inclinations. Choosing an homogeneous sampling of the initial conditions also avoids giving too much importance to regions protected by MMR owing to particular combinations of angles while another choice of angles would have given an unstable orbit.

4.1 Numerical setup

We ran numerical simulations using the symplectic scheme ABAH1064 from Farrés et al. (2013). We chose our initial conditions such that

• the outerplanet semimajor axis a₂ is fixed at 1 au;

• the AMD and inner planet semimajor axis a₁ are chosen such that, we have a regular grid in the plane $(\alpha ,\sqrt{\mathcal{C}})$. Such a scaling in $\mathcal{C}$ is chosen to have an approximately uniform distribution in terms of eccentricities and inclination;

• the AMD is on average equipartitioned between the eccentricity and inclination degrees of freedom;

• the inclinations are chosen such that the angular momentum is on the z-axis;

• the angles are chosen randomly;

• the star mass is taken as 1 M_⊙ and the planets masses do not vary for each grid of initial conditions.

We then integrate each initial condition for 500 kyr using a time-step of 10⁻³ yr. The numerical integration is stopped if the planets approach each other by less than a quarter of their mutual Hill radius, if a planet reaches 10⁻² or 20 AU or if the relative variation of energy is higher than 10⁻⁸.

In order to measure the stability of a system, we use frequency map analysis (Laskar 1990, 1993). Our criterion is based on the relative variation of the main frequencies in the quasiperiodic best fit. More precisely, let $n_k^{(i)}$ (respectively $n_k^{(f)}$) be the frequency obtained by frequency analysis for the planet k for the first (resp. last) 100 kyr of integration. We consider an orbit to be chaotic if $$\delta n={\displaystyle \underset{k}{\max }}\left|\frac{n_k^{(f)}-n_k^{(i)}}{n_k^{(i)}}\right|$$(31)

is greater than 10⁻⁴. The chosen threshold is such that the variation of semimajor axis changes by about 1% in a few Gyr if we assume a constant diffusion process as it would happen for a random walk. Indeed, δn measures the variation of frequency over 500 kyr. If the diffusion rate remains constant, a variation on the order of 1% on average needs a time 10 000 times larger, i.e., 5 Gyr.

If the integration time is shorter due to a collision or ejection, we set δn to 1. Since we randomly draw most of the initial parameters, we bin the results into a two-dimensional grid in $(\alpha ,\sqrt{\mathcal{C}})$ and average the frequency variation in each bin.

4.2 Results

We first integrate 1 00 000 initial conditions on a uniform grid with α taking values from 0.5 to 1 and $\mathcal{C}$ from 0 to 0.1. The masses of the two planets are equal to 0.5 × 10⁻⁵ M_⊙, such that ε = 10⁻⁵ and γ = 1. In this simulation, 78.4% of the orbits survive up to 500 kyr, 21.1% end up in a collision between the two planets, and 0.5% ofthe integrations are stopped because of the nonconservation of energy due to an unresolved close encounter.

The results of the frequency analysis are shown in Fig. 4. We see that the chaotic region is well constrained by the Hill curve $C_c^{\text{H}}$. Indeed, very few orbits with $\mathcal{C}<C_c^{\text{H}}$ appear to be chaotic. The region where Hill stable orbits ($\mathcal{C}<C_c^{\text{H}}$) are chaotic seems restricted to the region, where $C_{c^{\text{MMR}}}<\mathcal{C}<C_c^{\text{H}}$ around α ≃ 0.94 and low $\mathcal{C}$, i.e., for orbits experiencing MMR overlap (a zoomed view of Fig. 4 is given in Fig. 5). The behavior of planets initially in this region was already discussed in (Deck et al. 2013). Orbits that are Hill unstable ($\mathcal{C}>C_c^{\text{H}}$) appear to be largely chaotic up to some resonant islands situated at α ≃1 (co-orbital resonance) and near the 3:2 and 4:3 resonances (α ≃ 0.76 and 0.82). We also see that for larger separations (α ≲ 0.6), orbits are less chaotic. However, it is probable that for longer integration times, these orbits would end up unstable.

In order to highlight how the Hill criterion separates the chaotic orbits from the stable orbits, we can plot the fraction of regular orbits as a function of $(p/a)/{(p/a)}_c$ as suggested in Barnes & Greenberg (2006). We see in Fig. 6 that there is a sharp limit at $p/a={(p/a)}_c$. All Hill stable integrations go to 500 kyr and very few are chaotic. For Hill unstable orbits ($\mathcal{C}>C_c^{\text{H}}$), we see a slight increase of regular and surviving orbits with p∕a but the decrease is not as significant as the change at the Hill limit. On average 64.9% of the Hill unstable orbits survive and 52.3% are regular.

We plot in Fig. 7 the average time of the integrations as a function of initial α and $\mathcal{C}$. The initial conditions are binned in a 160 × 75 grid. We see in this figure that the average lifetime of a system is almost always smaller than 500 kyr in the Hill unstable region. We also remark that the average lifetime is less than 50 kyr for a system with α ≳0.9 and not too large $\mathcal{C}$ but increases for wider separations or greater AMD. Indeed, for higher $\mathcal{C}$, the planets are on average initially further away from crossing orbits because the initial choice of eccentricities and inclinations is random.

Fig. 4 Frequency variation δn (31) for 1 00 000 initial conditions binned in a 160 × 75 grid with ε = 10−5 and γ = 1. The black curve is the Hill critical AMD $C_c^{\text{H}}$ (30), the purple curve is the collisional critical AMD Cc (Laskar & Petit 2017), and the green curve is the critical AMD obtained from the overlap of first-order MMR $C_{c^{\text{MMR}}}$ (Petit et al. 2017). Each bin is average over about eight initial conditions.

Fig. 5 Zoomed view of Fig. 4 for 0.9 < α < 1 and $\sqrt{\mathcal{C}}<0.07$. In the regionwhere $C_{c^{\text{MMR}}}<\mathcal{C}<C_c^{\text{H}}$, we see that orbits are chaotic but Hill stable.

Fig. 6 Proportion of orbits that survived for 500 kyr (green curve) and nonchaotic orbits (δn < 10−4; red curve) asa function of $(p/a)/{(p/a)}_c$.$(p/a)/{(p/a)}_c>1$ means that the initial condition verifies the Hill criterion. Above the main figure, we added azoomed view of the upper part of the plot. The light blue histogram represents the number of initial conditions in each of the 300 bins used to compute the fractions plotted.

Fig. 7 Average lifetime of the system as a function of α and $\mathcal{C}$. The parameters are similar to Fig. 4. Each bin represents an average over eight initial conditions. Dark blue implies that all integration ended at 500 kyr.

4.3 Influence of the masses

To highlight the role of the mass of the planets, we ran another simulation with the same method but with ε = 10⁻³ and γ = 1. We integrate 40 000 initial conditions and carry out a similar analysis. After 500 kyr, 54.6% of the systems leads to a collision between the two planets, 2% to an ejection, and 0.5% is stopped due to nonconservation of energy. If we consider only the Hill unstable systems ($\mathcal{C}>C_c^{\text{H}}$), only 24.0% have survived and 15.7% are regular. As expected, larger masses lead to much more unstable systems.

The results of the frequency analysis are plotted in Fig. 8. We observe that the 3:2 MMR is still very stable in comparison to the surrounding regions. The 2:1 MMR appears to create an area of moderate chaos in the Hill stable region but is not as marked in the Hill unstable part. The co-orbital resonance also appears more unstable.

We see in Fig. 9 that some of the Hill stable systems do not survive for 500 kyr. These systems have led to an ejection of the outer planet. As a result, the Hill stability line appears to be more porous in comparison to the case ε = 10⁻⁵.

We also see that a larger percentage of Hill stable orbits appear to be chaotic. This can be quantified, thanks to Fig. 10. We see that the proportion of regular orbits drops before it reaches the Hill stability limit. However, the more important change of behavior is still at the Hill stability limit.

Fig. 8 Frequency variation δn (31) for 40 000 initial conditions binned in a 100 × 50 grid with ε = 10−3 and γ = 1. The three curves are similar to Fig. 4.

Fig. 9 Average lifetime of the system as a function of α and $\mathcal{C}$. The parameters are similar to Fig. 8. We see that some Hill stable initial conditions stopped before 500 kyr owing to ejections.

Fig. 10 Proportion of orbits that survived for 500 kyr (green curve) and nonchaotic orbits (δn < 10−4; red curve) asa function of $(p/a)/{(p/a)}_c$.$(p/a)/{(p/a)}_c>1$ means that the initial condition verifies the Hill criterion. The light blue histogram represents the number of initial conditions in each of the 150 bins used to compute the fractions plotted.

5 Conclusions

In a two-planet system, Hill stability is a topological limit that forbids close encounters between the outer planet and the inner planet or the star. If verified, the system remains stable if the outer planet does not escape or if the inner planet does not collide with the star. Moreover, since a minimal distance between planets is imposed, the planet perturbation remains moderate and the system is most likely regular.

We have generalized Gladman’s Hill stability criterion and have shown that it is natural to express the Hill criterion in the AMD framework. Indeed, we obtained a simple expression for this criterion (Proposition 3) with only an expansion in ε (Eq. (1)). Moreover, it is easy to recover all former published criteria as particular cases of this expression. Because of its formulation as a function of the AMD, the expression (26) is valid in the general spatial case for any value of eccentricities and inclinations. Moreover, our Hill criterion is accurate even for very different planet masses. We also highlight that the AMD and the semimajor axis ratio α are the main parameters to consider in a stability study.

We show that the Hill stability allows us to give an accurate stability limit up to large orbital separations. The sharp change of behavior at the Hill stability limit has already been studied in Barnes & Greenberg (2006) or in Deck et al. (2013). Nevertheless, our numerical integrations confirm it for a much larger range of α and AMD with randomized initial conditions for the parameters not taken into account.

Our simulations for large planet masses also show that the expansion in ε is valid even for larger planets and that the Hill stability accurately segregates between regular and short-lived initial conditions. However, it appears that for large mass values, ejections cannot be neglected and a model should be developed to understand this further behavior.

As shown in several works on tightly packed systems (Chambers et al. 1996; Pu & Wu 2015), Gladman’s Hill criterion is no longer adapted in these cases. Such a sharp limit between almost eternal and short-lived systems no longer exists. Instead, it appears that there exists a scaling between the initial orbital separation and the time of instability, wider separated orbits becoming unstable after a longer time. In the cited works, the empirical stability criteria give stability spacing as a function of the mutual radius $$\frac{a_{n+1}-a_n}{a_{n+1}}>K{R}_{\text{H}},$$(32)

where $R_{\text{H}}={(\epsilon /3)}^{1/3}$ is the Hill radius. In the context of multiplanetary systems, an analytical work on long-term stability is still necessary. This will be the subject of future work.

Acknowledgements

This work was granted access to the HPC resources of MesoPSL financed by the Region Ile de France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d’Avenir supervised by the Agence Nationale pour la Recherche. We acknowledge the support of the PNP and the CS of the Observatoire de Paris.

Appendix A Computation of R at the Lagrange points

The function $R(x,y)={(\rho /\nu )}^2$ defined in Eq. (16) admits three saddle points situated on the x-axis, which are the Lagrange points L₁, L₂, and L₃. If x_j is the abscissa of the point L_j, we have 0 < x₁ < 1, x₂ > 1, and x₃ < 0. The x_j quantities only depends on the mass ratio ε and on $$\zeta =\frac{\gamma }{\gamma +1}=\frac{m_1}{m_1+m_2}.$$(A.1)

They are the roots of the three different polynomial equations: $$\begin{array}{l}(1+\zeta \epsilon )x_1^5-(2+3\zeta \epsilon )x_1^4\\ +(1+2\epsilon +\zeta \epsilon )x_1^3-(1+3\epsilon (1-\zeta ))x_1^2\\ +(2+3\epsilon (1-\zeta ))x_1-(1+\epsilon (1-\zeta ))=0,\\ (1+\zeta \epsilon )x_2^5-(2+3\zeta \epsilon )x_2^4\\ +(1+3\zeta \epsilon )x_2^3-(1+3\epsilon -\zeta \epsilon ))x_2^2\\ +(2+3\epsilon (1-\zeta ))x_2-(1+\epsilon (1-\zeta ))=0,\\ (1+\zeta \epsilon )x_3^5-(2+3\zeta \epsilon )x_3^4\\ +(1+3\zeta \epsilon )x_3^3+(1+3\epsilon (1-\zeta ))x_3^2\\ +(2+3\epsilon (1-\zeta ))x_3+(1+\epsilon (1-\zeta ))=0,\end{array}$$

At up to terms of order ε^4∕3, we have $$\begin{array}{ll}{x}_1\hfill & =1-{\left(\frac{\epsilon }{3}\right)}^{1/3}+\frac{\gamma +2}{3(\gamma +1)}{\left(\frac{\epsilon }{3}\right)}^{2/3}+\frac{(\gamma -2)\epsilon }{27(\gamma +1)}+\text{O}({\epsilon }^{4/3}),\hfill \\ x_2\hfill & =1-{\left(\frac{\epsilon }{3}\right)}^{1/3}+\frac{\gamma +2}{3(\gamma +1)}{\left(\frac{\epsilon }{3}\right)}^{2/3}-\frac{(\gamma -2)\epsilon }{27(\gamma +1)}+\text{O}({\epsilon }^{4/3}),\hfill \\ x_3\hfill & =-1+\frac{7}{12}\frac{\gamma -1}{\gamma +1}\epsilon +\text{O}({\epsilon }^2).\hfill \end{array}$$(A.5)

We can then evaluate R (Eq. (16)) at those points and we have at order, ε^4∕3: $$\begin{array}{ll}R(L_1)\hfill & =1+\frac{3^{4/3}{\epsilon }^{2/3}\gamma }{{(\gamma +1)}^2}-\frac{(11+7\gamma )\gamma \epsilon }{3{(\gamma +1)}^3}+\text{O}({\epsilon }^{4/3}),\hfill \\ R(L_2)\hfill & =1+\frac{3^{4/3}{\epsilon }^{2/3}\gamma }{{(\gamma +1)}^2}-\frac{(11\gamma +7)\gamma \epsilon }{3{(\gamma +1)}^3}+\text{O}({\epsilon }^{4/3}),\hfill \\ R(L_3)\hfill & =1+\frac{2\epsilon \gamma }{{(\gamma +1)}^2}+\text{O}({\epsilon }^2).\hfill \end{array}$$(A.6)

Appendix B Expansion of $C_c^{\text{Ex}}$ and h₁

In Sect. 2, the Hill stability criterion (26) is obtained by the expansion at the leading order in ε of $$F={(\gamma +1)}^{3/2}\sqrt{\frac{R(L_1){(1+\epsilon \gamma /{(\gamma +1)}^2)}^3}{(1+\epsilon )\left(\gamma /\alpha +1+h_1\right)}}.$$(B.1)

In F, the main term depending on ε comes from the expansion of $$R(L_1)=1+\frac{3^{4/3}\gamma {\epsilon }^{2/3}}{{(\gamma +1)}^2}-\frac{(11+7\gamma )}{3(\gamma +1)}\frac{\gamma \epsilon }{{(\gamma +1)}^2}+\text{O}({\epsilon }^{4/3}).$$(B.2)

However, wealso need to make sure that h₁ remains small in comparison to this term. Thus, the expansion of F requires an estimate of h₁ with respect to ε and γ. As explained in Sect. 2, h₁ is the renormalized perturbation part of the Hamiltonian (3): $$h_1=-\frac{2{\Lambda }_2^2}{m_2^3{\mu }^2}\left(\frac{1}{2}\frac{{\Vert {\stackrel{}{r\sim }}_1+{\stackrel{}{r\sim }}_2\Vert }^2}{m_0}-\frac{\mathcal{G}{m}_1{m}_2}{r_{12}}\right).$$(B.3)

If we note ${\stackrel{.}{r}}_j={\stackrel{}{r\sim }}_j/m_j$ and simplify the expression (B.3), we obtain $h_1=h_1^{\text{T}}+h_1^{\text{P}}$, where $$h_1^{\text{T}}=-\frac{\epsilon a_2}{\mu }\frac{\Vert \gamma {\stackrel{.}{r}}_1+{\stackrel{.}{r}}_2{\Vert }^2}{\gamma +1},\text{\hspace{1em}and\hspace{1em}}{h}_1^{\text{P}}=\frac{2a_2}{r_{12}}\frac{\epsilon \gamma }{\gamma +1}.$$(B.4)

As one cansee in Fig. 2, the value of ε is only relevant for small values of the critical AMD $C_c^{\text{Ex}}$, i.e., for close planets. Using the expansion in 1 − α (27), we can estimatethat we need to compute h₁ for systems such that $\mathcal{C}$ is of order ε^2∕3. This corresponds to systems with eccentricities of order ε^1∕3. We can therefore use the circular approximation in our estimation of h₁. In particular, we have r_j = a_j(1 + O(ε^1∕3)).

We first consider the term coming from the gravitational interaction between the two planets, $h_1^{\text{P}}$. If the system is Hill stable, the distance r₁₂ is greater than the radius ofthe Hill sphere $S_{{\text{H}}_{\text{1}}}$: $$\begin{array}{ll}{r}_{12}>r_1{\max }_{j=1,2}|1-x_j|\hfill & =r_1{(\epsilon /3)}^{1/3}+\text{O}({\epsilon }^{2/3})\hfill \\ \hfill & =a_1{(\epsilon /3)}^{1/3}+\text{O}({\epsilon }^{2/3}),\hfill \end{array}$$(B.5)

where x_j is defined in (A.5). For all times, we therefore have $$h_1^{\text{P}}=\frac{2a_2}{r_{12}}\frac{\epsilon \gamma }{\gamma +1}\le \frac{2\times 3^{1/3}}{\alpha }{\epsilon }^{2/3}\frac{\gamma }{\gamma +1}+\text{O}(\epsilon ).$$(B.6)

The gravitational potential term $h_1^{\text{P}}$ is at most of the same order as the leading term in ε of R(L₁). However, we can always choose to estimate the energy and actions values when the two planets are far from each other. In this case, a₂ ∕r₁₂ = O(1) and $h_1^{\text{P}}$ is linear in ε. From now on, we assume that we have $$h_1^{\text{P}}=\frac{\epsilon \gamma }{\gamma +1}{p}_{12},$$(B.7)

with p₁₂ = O(1).

Moreover, F is a decreasing function of h₁, so $C_c^{\text{Ex}}$ is an increasing function of h₁. Since $h_1^{\text{P}}$ is positive, neglecting it is equivalent to have a more conservative criterion.

Let us now consider the kinetic term $h_1^{\text{T}}$. We develop (B.4) for almost circular orbits. In this limit, we have $\Vert {\stackrel{.}{r}}_j\Vert =\sqrt{\mu /a_j}+\text{O}({\epsilon }^{1/3})$. We have $$h_1^{\text{T}}=-\frac{\epsilon }{\gamma +1}\left(\frac{{\gamma }^2}{\alpha }+1+\frac{2\gamma }{\sqrt{\alpha }}\cos ({\lambda }_1-{\lambda }_2)\right)+\text{O}({\epsilon }^{4/3}).$$(B.8)

Therefore, $h_1^{\text{T}}$ is always linear in ε. Combining the estimations (B.7) and (B.8), we see that h₁ = O(ε).

We can now use the expansions of h₁, (B.8) and (B.7), to obtain the expansion of F in function of ε and γ. Because of the term $3^{4/3}{\epsilon }^{2/3}\gamma /{(\gamma +1)}^2$ from R(L₁), we do not keep terms of order $\text{O}(\epsilon \gamma /{(\gamma +1)}^2)$. We keep all other terms depending on ε up to the order O(ε^4∕3). We obtain $$F=\sqrt{\frac{\alpha {(\gamma +1)}^3}{\gamma +\alpha }\left(1+\frac{3^{4/3}{\epsilon }^{2/3}\gamma }{{(\gamma +1)}^2}\right)D}+\text{O}\left(\frac{\epsilon \gamma }{{(\gamma +1)}^2},{\epsilon }^{4/3}\right),$$(B.9)

where $D=(1-\epsilon )\left(1-\frac{\alpha h_1}{\gamma +\alpha }\right)$ and $$\text{O}\left(\frac{\epsilon \gamma }{{(\gamma +1)}^2},{\epsilon }^{4/3}\right)=\text{O}\left(\frac{\epsilon \gamma }{{(\gamma +1)}^2}\right)+\text{O}\left({\epsilon }^{4/3}\right).$$

Let us develop D on the same order than F. We have $$\begin{array}{ll}D\hfill & =1-\epsilon -\frac{\alpha h_1^{\text{T}}}{\gamma +\alpha }+\text{O}\left(\frac{\epsilon \gamma }{{(\gamma +1)}^2},{\epsilon }^{4/3}\right)\hfill \\ \hfill & =1+\epsilon \left(\frac{{\gamma }^2+\alpha }{(\gamma +1)(\gamma +\alpha )}-1\right)+\text{O}\left(\frac{\epsilon \gamma }{{(\gamma +1)}^2},{\epsilon }^{4/3}\right)\hfill \\ \hfill & =1-\frac{\epsilon \gamma }{{(\gamma +1)}^2}\frac{(\gamma +1)(\alpha +1)}{\gamma +\alpha }+\text{O}\left(\frac{\epsilon \gamma }{{(\gamma +1)}^2},{\epsilon }^{4/3}\right)\hfill \\ \hfill & =1+\text{O}\left(\frac{\epsilon \gamma }{{(\gamma +1)}^2},{\epsilon }^{4/3}\right).\hfill \end{array}$$(B.10)

We see that D has no contribution to F at the considered orders. Therefore, we can write $$F=\sqrt{\frac{\alpha {(\gamma +1)}^3}{\gamma +\alpha }\left(1+3^{4/3}\frac{{\epsilon }^{2/3}\gamma }{{(\gamma +1)}^2}\right)}+\text{O}\left(\frac{\epsilon \gamma }{{(\gamma +1)}^2},{\epsilon }^{4/3}\right).$$(B.11)

As an immediate consequence of (B.11), the expression proposed in Eq. (26) remains valid in the case of a very uneven mass distribution between the two planets (e.g., for O (ε^1∕3) < γ < O(ε^−1∕3)).

References

All Figures

Fig. 1 Some levels of the function R defined in (16) for ε = 10−3 and γ = 1 in the $\mathcal{P}$ plane. The two red points correspond to the Lagrange points L1 and L2 and the three orange points to L3, L4, and L5. The orange-filled area corresponds to the region where R(x, y) < R(L1). The star S is at the origin, the planet P1 at (1,0), and (x,y) are the coordinates of the second planet. The green region $S_{{\text{H}}_{\text{1}}}$represents the Hill sphere of the first planet.

In the text

Fig. 2 Comparison of the right-hand sides of the inequalities – Eq. (25) (orange), Eq. (26) (green), and Eq. (27) (red) – as a function of α for γ = 1 and ε = 10−5 (upper three curves) or ε = 10−3 (lower three curves). The black curve is the critical AMD for the collision condition (Laskar & Petit 2017). The green and orange curves are on top of each other (see Fig. 3). The critical AMD from MMR overlap$C_{c^{\text{MMR}}}$ and circular criterion of Gladman (1993) are also plotted for comparison.

In the text

	Fig. 3 Maximum difference between $C_c^{\text{Ex}}$ (Eq. (25)) and $C_c^{\text{H}}$ (Eq. (30)) on the right-hand sides of Eqs. (25) and (26), respectively, as a function of ε for various values of γ = m1∕m2. In Eq. (25), the term h1 is evaluated using the approximation of the kinetic term given in Appendix B.
In the text

Fig. 4 Frequency variation δn (31) for 1 00 000 initial conditions binned in a 160 × 75 grid with ε = 10−5 and γ = 1. The black curve is the Hill critical AMD $C_c^{\text{H}}$ (30), the purple curve is the collisional critical AMD Cc (Laskar & Petit 2017), and the green curve is the critical AMD obtained from the overlap of first-order MMR $C_{c^{\text{MMR}}}$ (Petit et al. 2017). Each bin is average over about eight initial conditions.

In the text

	Fig. 5 Zoomed view of Fig. 4 for 0.9 < α < 1 and $\sqrt{\mathcal{C}}<0.07$. In the regionwhere $C_{c^{\text{MMR}}}<\mathcal{C}<C_c^{\text{H}}$, we see that orbits are chaotic but Hill stable.
In the text

Fig. 6 Proportion of orbits that survived for 500 kyr (green curve) and nonchaotic orbits (δn < 10−4; red curve) asa function of $(p/a)/{(p/a)}_c$.$(p/a)/{(p/a)}_c>1$ means that the initial condition verifies the Hill criterion. Above the main figure, we added azoomed view of the upper part of the plot. The light blue histogram represents the number of initial conditions in each of the 300 bins used to compute the fractions plotted.

In the text

	Fig. 7 Average lifetime of the system as a function of α and $\mathcal{C}$. The parameters are similar to Fig. 4. Each bin represents an average over eight initial conditions. Dark blue implies that all integration ended at 500 kyr.
In the text

	Fig. 8 Frequency variation δn (31) for 40 000 initial conditions binned in a 100 × 50 grid with ε = 10−3 and γ = 1. The three curves are similar to Fig. 4.
In the text

	Fig. 9 Average lifetime of the system as a function of α and $\mathcal{C}$. The parameters are similar to Fig. 8. We see that some Hill stable initial conditions stopped before 500 kyr owing to ejections.
In the text

	Fig. 10 Proportion of orbits that survived for 500 kyr (green curve) and nonchaotic orbits (δn < 10−4; red curve) asa function of $(p/a)/{(p/a)}_c$.$(p/a)/{(p/a)}_c>1$ means that the initial condition verifies the Hill criterion. The light blue histogram represents the number of initial conditions in each of the 150 bins used to compute the fractions plotted.
In the text