Antoine C Petit
email: antoine.petit@obspm.fr

Jacques Laskar

Gwenaël Boué

Hill stability in the AMD framework

Keywords: Celestial mechanics -Planets and satellites, general -Planets and satellites, dynamical evolution and stability

 inequality, a topological boundary can forbid close encounters between the two planets for infinite time. A system is said Hill stable if it verifies this topological condition. Hill stability is widely used in the study of extra solar planets dynamics. However people often use the coplanar and circular orbits approximation.

In this paper, we explain how the Hill stability can be understood in the framework of Angular Momentum Deficit (AMD). In the secular approximation, the AMD allows to discriminate between a priori stable systems and systems for which a more in depth dynamical analysis is required. We show that the general Hill stability criterion can be expressed as a function of only the semi major axes, the masses and the total AMD of the system. The proposed criterion is only expanded in the planets-to-star mass ratio ε and not in the semi-major axis ratio, in eccentricities nor in the mutual inclination. Moreover the expansion in ε remains excellent up to values of about 10 -3 even for two planets with very different mass values. We performed numerical simulations in order to highlight the sharp change of behaviour between Hill stable and Hill unstable systems. We show that Hill stable systems tend to be very regular whereas Hill unstable ones often lead to rapid planet collisions. We also remind that Hill stability does not protect from the ejection of the outer planet.

Introduction

The chaotic nature of planetary dynamics has made necessary the development of stability criteria 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 exoplanets dynamics, the orbital parameters are often not known with great precision, making impossible to conduct a precise dynamical study. For multi-planetary 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) has 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 1997Laskar , 2000)). The AMD is a weighted sum of the planets' eccentricities and mutual inclinations and can be interpreted as a dynamical temperature of the planetary system. The AMD is conserved at all order 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, the AMD-stability. AMD-stable systems 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 in order 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 since, have 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 show in (Petit et al. 2017) that this MMR overlap criterion can be expressed solely as a function of the AMD, semi-major 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 the 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, the Hill stability. Based on a work by Sundman (1912) on the moment of inertia in the three-body problem, Marchal & Saari (1975) note the existence of forbidden zones in the config-uration space. In (Marchal & Bozis 1982), they extend the notion of Hill stability to the general three-body problem and show 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;Marchal & Bozis 1982) have many applications outside of the Hill stability. One can 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 triple close approach for bounded orbits (Laskar & Marchal 1984). Marchal & Bozis (1982) present 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's paper also proposes some criteria for eccentric orbits in some particular configurations of masses. Gladman's result was refined in order to cover other situations as, for example, the case of inclined orbits (Veras & Armitage 2004). See (Georgakarakos 2008) for 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) analyze the differences between the MMR overlap and Hill stability criteria. They remark that there exists an area where orbits are chaotic due to the overlap of MMR but long lived due to Hill stability.

In this paper, we show that Marchal & Bozis (1982) Hill stability criterion fits extremely well in the AMD-stability framework. In section 2, we derive a criterion for Hill stability solely expressed as a function of the total AMD, the semi-major axes and the 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;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 semi-major axis ratio.

Hill stability in the 3-body problem

Generalised Hill curves

Let us use the formalism proposed in (Marchal & Saari 1975;Marchal & Bozis 1982) for the definition of the Hill regions for the general 3-body problem. We mainly consider the planetary case and therefore adapt slightly their notations to this particular problem. Let us consider two planets of masses m 1 , m 2 orbiting a star of mass m 0 . Let G be the constant of gravitation, µ = Gm 0 ,

ε = m 1 + m 2 m 0 , (1)
the planet mass to star mass ratio and

γ = m 1 m 2 . (2)
As in Laskar & Petit (2017), we use the heliocentric canonical coordinates (r j , r j ) j=1,2 . In those coordinates, the Hamiltonian is

H = 2 j=1         1 2 r j 2 m j - µm j r j         + 1 2 r1 + r2 2 m 0 - Gm 1 m 2 r 12 , (3)
where r 12 = r 1r 2 . We denote G the total angular momentum of the system

G = 2 j=1 r j × r j , (4)
and G its norm. From now, we assume G to be aligned along the z-axis. We also define

Λ j = m j √ µa j , (5)
where a j is the semi-major 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(Laskar , 2000) ) has for expression

C = 2 j=1 Λ j 1 -1 -e 2 j cos(i j ) . (6)
Following (Marchal & Bozis 1982), we define a generalised semi-major axis

a = - GM * 2H , (7)
and a generalised semi-latus rectum

p = m 0 (1 + ε) GM * 2 G 2 , (8)
where

M * = m 0 m 1 + m 0 m 2 + m 1 m 2 = m 2 0 ε 1 + ε γ (1 + γ) 2 . ( 9
)
a and p are the two length units that can be built from the first integrals H and G.

We finally define two variable lengths. First, ρ the mean quadratic distance

M * ρ 2 = m 0 m 1 r 2 1 + m 0 m 2 r 2 2 + m 1 m 2 r 2 12 , (10)
that is proportional to the moment of inertia I computed in the centre of mass frame (Marchal & Bozis 1982)

I = 1 2 M * ρ 2 m 0 + m 1 + m 2 . ( 11
)
We also define ν, the mean harmonic distance

M * ν = m 0 m 1 r 1 + m 0 m 2 r 2 + m 1 m 2 r 12 , (12)
which is proportional to the potential energy

U = - GM * ν . ( 13
)
For a system with given H and G, some configurations of the planets are forbidden. Indeed, the value of the ratio ρ/ν is constrained by the inequality (Marchal & Saari 1975)

ρ ν ≥ ρ 2a + p 2ρ , (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 ρ = √ ap. We therefore have the inequality (Marchal & Bozis 1982)

ρ 2 ν 2 ≥ p a = - 2m 0 (1 + ε) G 2 M * 3 HG 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 (ρ/ν) 2 . Since this ratio only depends on the ratios of mutual distances, we can always place ourselves in the plane generated by the three bodies. We can also choose to place the first planet on the x-axis and normalize the lengths by r 1 . Let us call this plane P and note (x, y) the coordinates of the second planet. In the plane P (see Figure 1), the star S is at the origin, the first planet P 1 is situated at the point (1, 0). For the planar restricted 3-body problem, this reduction is equivalent to study the dynamics in the corotating frame. We note

R(x, y) = ρ ν 2 . ( 16
)
The shape of the function R in P only depends on the mass distribution, e.g., the two ratios ε and γ. In Figure 1, we can see level curves of the function R plotted in P for ε = 10 -3 and γ = 1. R is minimal and equal to 1 at the two Lagrange points L 4 and L 5 and goes to +∞ for r 2 → 0, r 2 -r 1 → 0 or r 2 → +∞. The function has 3 saddle nodes at the Lagrange points L 1 , L 2 and L 3 . We give the method to compute the position of the Lagrange points in appendix A. Let us denote respectively x 1 , x 2 and x 3 their abscissa. At the lowest order in ε we have

x 1 = 1 - ε 3 1/3 + O(ε 2/3 ), x 2 = 1 + ε 3 1/3 + O(ε 2/3 ), (17)
x 3 = -1 + 7 12 γ -1 γ + 1 ε + O(ε 2 ).
At the lowest order in ε, the value of R at the Lagrange points are

R(L 1 ) = 1 + 3 4/3 ε 2/3 γ/(γ + 1) 2 + O(ε) R(L 2 ) = 1 + 3 4/3 ε 2/3 γ/(γ + 1) 2 + O(ε), (18)
R(L 3 ) = 1 + 2εγ/(γ + 1) 2 + O(ε 2 ).
R(L 1 ) and R(L 2 ) have the same first order term but differ in the expansion of higher order. Indeed, if

m 0 ≥ m 1 ≥ m 2 , we have R(L 1 ) ≥ R(L 2 ) (Marchal & Bozis 1982). From now, let us assume R(L 1 ) ≥ R(L 2 ) (if not, we can just substitute R(L 2 ) to R(L 1 ) in further equations). x y L 1 L 2 L 3 L 4 L 5 S P 1 S H1
Fig. 1. We represent in the P plane, some levels of the function R defined in ( 16) for ε = 10 -3 and γ = 1. The two red points correspond to the Lagrange points L 1 and L 2 and the three orange points to L 3 , L 4 and L 5 . The orange-filled area corresponds to the region where R(x, y) < R(L 1 ). The star S is at the origin, the planet P 1 at (1,0) and (x,y) are the coordinates of the second planet. The green region S H 1 represents the Hill sphere of the first planet.

For p/a ≥ R(L 1 ), the accessible domain is split in three parts, the Hill sphere of the star S which is around the origin, the first planet Hill sphere1 S H 1 (in green on Figure 1) and the outer region. In this case, if the second planet is not initially inside S H 1 , it will never be able to enter this region. However, if P 2 is in the outer region, the Hill stability cannot constrain the possibility of ejection. Similarly, if P 2 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 non-collision criterion for an infinite time. Marchal and Bozis called it the Hill stability.

Proposition 1 (Marchal & Bozis 1982). Let us consider a negative energy three-body problem with a body S of mass m 0 and two others P 1 and P 2 of mass m 1 and m 2 such that m 0 ≥ m 1 ≥ m 2 . We place ourselves in the P plane defined by S , P 1 and P 2 (Figure 1). If P 2 is not initially inside the Hill sphere S H 1 of P 1 , the system is Hill stable if

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.

AMD condition for Hill stability

The total energy of the system can be written

H = - m 3 2 µ 2 2Λ 2 2 γ α + 1 + h 1 , (20)
where α = a 1 /a 2 and

h 1 = - 2Λ 2 2 m 3 2 µ 2 1 2 r1 + r2 2 m 0 - Gm 1 m 2 r 12 . ( 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 = Λ 2 γ √ α + 1 -C , (22)
where

C = C Λ 2 = γ √ α 1 -1 -e 2 1 cos i 1 + 1 -1 -e 2 2 cos i 2 (23)
is the relative AMD (Laskar & Petit 2017). Combining ( 15), ( 20) and ( 22), we obtain

p a = (1 + ε) (γ + 1) 3 (1 + εγ/(γ + 1) 2 ) 3 γ α + 1 + h 1 γ √ α + 1 -C 2 . ( 24
)
The Hill stability criterion (19) can be rewritten without any approximation as a condition on C and we have the following formulation of the Hill stability.

Proposition 2 (Hill stability). With the hypotheses of the proposition 1, assuming the elliptical elements can be defined ( i.e. both Keplerian energies are negative), a system is Hill stable if

C < C Ex c = γ √ α+1-(γ+1) 3/2 R(L 1 )(1 + εγ/(γ + 1) 2 ) 3 (1 + ε) (γ/α + 1 + h 1 ) , ( 25
)
where C is the relative AMD ( 23) and h 1 the normalized perturbation part (21).

The inequality ( 25) is equivalent to the proposition 1 but we isolated on the left hand side the contribution from the AMD. Up to the perturbation term h 1 , the right hand side of (25) only depends on the masses and the semi-major 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 γ.

Proposition 3 (Hill stability, planetary case). For small enough ε, a two-planet system is Hill stable if the relative AMD C verifies the inequality

C < γ √ α + 1 -(1 + γ) 3/2 α γ + α 1 + 3 4/3 ε 2/3 γ (1 + γ) 2 + O(ε). (26)
As explained in appendix B, h 1 is of smaller order in ε and can be neglected if the criterion is verified. We want to stress out that the expression ( 26) is obtained with only an expansion in ε and does only depends of α, C and the masses of the bodies. The term of order ε 2/3 in (eq. 26) also depends on γ/(γ + 1) 2 , but we show in appendix B that ( 26) is still valid for γ 1 or γ 1.

Close planets approximation

Assuming 1 -α 1, and a small AMD value, further approximations can be made. At leading order in C , 1 -α and ε, the inequality (26) becomes

C < 3γ 8(γ + 1) (1 -α) 2 - 3 4/3 γ 2(γ + 1) ε 2/3 . ( 27
)
One can isolate 1 -α in this expression to obtain an approximate minimum spacing for Hill stable systems.

Proposition 4 (Hill stability, close planets case). For a close planets system, the minimum spacing criterion for Hill stability is

a 2 -a 1 a 2 = 1 -α > 4 × 3 1/3 ε 2/3 + 8 3 γ + 1 γ C . (28)
Gladman's eccentric criteria can be recovered from (28) if the AMD is developed under the assumptions (same planet masses, small or large and equal eccentricities) made in (Gladman 1993). However ( 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 -α > 2 × 3 1/6 ε 1/3 = 2.40ε 1/3 .
(29)

Comparison of the Hill criteria

We can compare the right-hand side of ( 25), ( 26), ( 27) and Gladman's circular approximation (29) to test how relevant are the approximations made here. In Figure 2, we plot the exact expression C Ex c (eq. 25, in green), the expansion in ε (eq. 26, in orange), the approximation for close planets (27, in red) as well as the minimum spacing for circular orbits (29, in blue). We see that the expansion in ε (26) cannot be distinguished from the exact curve (25). In order to better quantify this, we plot in Figure 3 the maximum difference between the two curves as a function of ε for various values of the mass ratio γ.

We see in figure 3 that for the range of ε used in planetary dynamics (typically from 10 -6 to 10 -3 ), the expression (26) developed in ε is accurate even for very uneven planet mass distribution. From now, we use (26) to define the Hill stability.

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 0.0 0.2 0.4 0.6 0.8 1.0 α -4 -3 -2 -1 0 log 10 C ε =10 -3 ε =10 -5 C c C MMR c C H c C Ex c
Developped in 1α Gladman circular spacing 25) and ( 26), as a function of ε for various values of γ = m 1 /m 2 . In (25), the term h 1 is evaluated using the approximation of the kinetic term given in appendix B obtained by finding the minimal AMD needed for a system to allow for collisions in the secular approximation, i.e., assuming the semi major axes are fixed. We completed in (Petit et al. 2017) the definition by excluding the configurations such that the overlap of first order MMR occurs. The modified critical AMD is hereafter noted

C MMR c .
We can use (26) to define a critical AMD C H c for the Hill stability

C H c = γ √ α + 1 -(1 + γ) 3/2 α γ + α 1 + 3 4/3 ε 2/3 γ (1 + γ) 2 . ( 30
)
A system will be Hill stable if its initial relative AMD is smaller than the initial critical AMD C H c . We see that for two planets, the Hill stability definition fits extremely well in the AMD-stability framework. We can also compare C H c to the previously proposed critical AMD. The collision critical AMD C c is plotted in Figure 2 with two values of C H c for ε = 10 -5 (resp. 10 -3 ). We can see that the Hill stability criterion is stricter (C H c < 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 C < C H c . As a result the two planets cannot approach each other by less than their mutual Hill radius for any variation of semi major 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 C < C c . Since we are not making any additional hypothesis on C , we have C H c < C c . The strict inequality comes from the positive minimal distance between the two planets.

As a comparison, we also plot in Figure 2

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) note 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 if 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 C and the ratio of semi-major 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 due to particular combinations of angles while another choice of angles would have given an unstable orbit.

Numerical set-up

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

the outer planet semi-major axis a 2 is fixed at 1 au, -the AMD and inner planet semi-major axis a 1 are chosen such that we have a regular grid in the plane (α, √ C ). Such a scaling in 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 1M 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 -3 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 -2 AU or 20 AU or if the relative variation of energy is higher than 10 -8 . In order to measure the stability of a system, we use the frequency map analysis (Laskar 1990(Laskar , 1993)). Our criterion is based on the relative variation of the main frequencies in the quasiperiodic best fit. More precisely, let be n (i) k (respectively n ( f ) k ) 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

δn = max k n ( f ) k -n (i) k n (i) k ( 31
)
is greater than 10 -4 . The chosen threshold is such that the variation of semi-major 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 of 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 draw randomly most of the initial parameters, we bin the results into a two-dimension grid in (α, √ C ) and average the frequency variation in each bin.

Results

We first integrate 100,000 initial conditions on a uniform grid with α taking values from 0.5 to 1 and C from 0 to 0.1. The masses of the two planets are equal to 0.5×10 -5 M , so that ε = 10 -5 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% of the integrations are stopped because of the non-conservation of energy due to an unresolved close encounter.

The 5 to see it more easily). The behavior of planets initially in this region was already discussed in (Deck et al. 2013). Orbits that are Hill unstable (C > C H c ) 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 ones, 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 Figure 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 (C > C H c ), 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 Figure 7 the average time of the integrations as a function of initial α and 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 C but increases for wider separations or greater AMD. Indeed for higher C , the planets are on average initially further away from crossing orbits because the initial choice of eccentricities and inclinations is random.

Influence of the masses

To highlight the role of the mass of the planets, we run another simulation with the same method but with ε = 10 -3 and γ = 1. We integrate 40,000 initial conditions and do a similar analysis. After 500 kyr, 54.6% of the systems lead to a collision between the two planets, 2% to an ejection and 0.5% are stopped due to non conservation of energy. If we consider only the Hill unstable systems (C > C H c ), 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 Figure 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 coorbital resonance also appears more unstable.

We see in Figure 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 -5 .

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

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 will remain stable if the outer planet does not escape or if the inner planet does not collide 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 generalised here Gladman's Hill stability criterion and have shown that it is natural to express the Hill criterion in the AMD framework. Indeed, we have 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 semi-major axis ratio α are the main parameters to consider in a stability study.

We show that the Hill stability allows to give an accurate stability limit up to large orbital separations. The sharp change of behaviour 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 randomised 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 behaviour.

As shown in several works on tightly packed systems (Chambers et al. 1996;Pu & Wu 2015), the Gladman's Hill criterion is not adapted anymore 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 a n+1a n a n+1 > KR H ,

where R H = (ε/3) 1/3 is the Hill radius.

In the context of multiplanetary systems, an analytical work on longterm stability is still necessary. This will be the subject of future work. As one can see on figure 2, the value of ε is only relevant for small values of the critical AMD C Ex c , i.e. for close planets. Using the expansion in 1 -α (27), we can estimate that we need to compute h 1 for systems such that C is of order ε 2/3 . It corresponds to systems with eccentricities of order ε 1/3 . We can therefore use the circular approximation in our estimation of h 1 . 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 P 1 . If the system is Hill stable, the distance r 12 is greater than the radius of the Hill sphere S H 1 r 12 > r 1 max j=1,2 |1x j | = r 1 (ε/3) 1/3 + O(ε 2/3 ) = a 1 (ε/3) 1/3 + O(ε 2/3 ), (B.5) where x j is defined in (A.5). For all times, we therefore have

h P 1 = 2a 2 r 12 εγ γ + 1 ≤ 2 × 3 1/3 α ε 2/3 γ γ + 1 + O(ε). (B.6)
The gravitational potential term h P 1 is at most of the same order as the leading term in ε of R(L 1 ). 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 2 /r 12 = O(1) and h P 1 is linear in ε. We will from now assume that we have Moreover, F is a decreasing function of h 1 , so C Ex c is an increasing function of h 1 . Since h P 1 is positive, neglecting it is equivalent to have a more conservative criterion.

Let us now consider the kinetic term h T 1 . We develop (B.4) for almost circular orbits. In this limit, we have ṙ j = µ/a j + O(ε 1/3 ). We have We can now use the expansions of h 1 , (B.8) and (B.7) to obtain the expansion of F in function of ε and γ. Because of the term 3 4/3 ε 2/3 γ/(γ + 1) 2 from R(L 1 ), we do not keep terms of order O(εγ/(γ + 1) 2 ). We keep all other terms depending on ε up to the order O(ε 4/3 ). We obtain As an immediate consequence of (B.11), the expression proposed in (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 )).

h T 1 = - ε γ + 1 γ 2 α + 1 + 2γ √ α cos(
F = α(γ + 1) 3 γ + α 1 + 3 4/3 ε 2/3 γ (γ +

Fig. 2 .Fig. 3 .

 23 Fig.2. Comparison of the right hand sides of the inequalities (25, orange), (26, green) and (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 eachother (see Figure3). The critical AMD from MMR overlap C MMR c andGladman's (1993) circular criterion are also plotted for comparison.

  the MMR critical AMD C MMR c . For small relative AMD C , C MMR c and C H c are almost identical.

  results of the frequency analysis are shown in Figure 4. We see that the chaotic region is well constrained by the Hill curve C H c . Indeed, very few orbits with C < C H c appear to be chaotic. The region where Hill stable orbits (C < C H c ) are chaotic seems restricted to the region where C MMR c < C < C H c around α 0.94 and low C i.e., for orbits experiencing MMR overlap (a zoom of Fig. 4 is given in Fig.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Frequency variation δn (31) for 100,000 initial conditions binned in a 160 × 75 grid with ε = 10 -5 and γ = 1. The black curve is the Hill critical AMD C H c (30), the purple curve is the collisional critical AMD C c (Laskar & Petit 2017) and the green one is the critical AMD obtained from the overlap of first order MMR C MMR c (Petit et al. 2017). Each bin is average over about 8 initial conditions.

Fig. 6 .

 6 Fig. 6. Proportion of orbits that survived for 500 kyr (green curve) and non chaotic orbits (δn < 10 -4 , red curve) as a function of (p/a)/(p/a) c . (p/a)/(p/a) c > 1 means that the initial condition verifies the Hill criterion. Above of the main figure, we added a zoom 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 .

 7 Fig. 7. Average lifetime of the system as a function of α and C . The parameters are similar to figure 4. Each bin represents an average over 8 initial conditions. A dark blue color implies that all integration ended at 500 kyr.

Fig. 8 .Fig. 9 .

 89 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 figure 4.

Fig. 10 .

 10 Fig. 10. Proportion of orbits that survived for 500 kyr (green curve) and non chaotic orbits (δn < 10 -4 , red curve) as a 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.

  we note ṙ j = r j /m j and simplify the expression (B.3), we obtain h 1 = h T 1 + h P 1

  λ 1 -λ 2 ) + O(ε 4/3 ). (B.8) Therefore, h T 1 is always linear in ε. Combining the estimations (B.7) and (B.8), we see that h 1 = O(ε).

D

  1) 2 D + O εγ (γ + 1) 2 , ε 4/3 , (B.9) where D = (1 -ε) 1 -αh 1 γ + α and O εγ (γ + 1) 2 , ε 4/3 = O εγ (γ + 1) 2 + O ε 4/3 .Let us develop D at the same order than F. We haveD = 1 -εαh T 1 γ + α + O εγ (γ + 1) 2 , ε 4/3 = 1 + ε γ 2 + α (γ + 1)(γ + α) -1 + O εγ (γ + 1) 2 , ε 4/3 Dhas no contribution to F at the considered orders. Therefore, we can writeF = α(γ + 1) 3 γ + α 1 + 3 4/3 ε 2/3 γ (γ + 1) 2 + O εγ (γ + 1) 2 , ε 4/3 .(B.11)

The Hill region is usually called the Hill sphere although it is not technically a sphere.

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 of the CS of the Observatoire de Paris

Appendix A: Computation of R at the Lagrange points

The function R(x, y) = (ρ/ν) 2 defined in (16) admits three saddle points situated on the x axis which are the Lagrange points L 1 , L 2 and L 3 . If x j is the abscissa of the point L j , we have 0 < x 1 < 1, x 2 > 1 and x 3 < 0. The x j quantities only depends on the mass ratio ε and on

They are the roots of the three different polynomial equations

At up to terms of order ε 4/3 , we have

+ O(ε 4/3 ), (A.5)

We can then evaluate R (eq. 16) at those points and we have at order,

Appendix B: Expansion of C Ex

c and h 1

In section 2, the Hill stability criterion ( 26) is obtained by the expansion at the leading order in ε of

In F, the main term depending on ε comes from the expansion of

However, we also need to make sure that h 1 remains small in comparison to this term. Thus, the expansion of F requires an estimate of h 1 with respect to ε and γ. As explained in section 2, h 1 is the renormalized perturbation part of the Hamiltonian (3)