EDP Sciences
Free Access
Issue
A&A
Volume 605, September 2017
Article Number A72
Number of page(s) 16
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201630022
Published online 12 September 2017

© ESO, 2017

1. Introduction

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

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

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

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

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

2. Angular momentum deficit

2.1. Planetary Hamiltonian

Let P0,P1,...,Pn be n + 1 bodies of masses m0,m1,...,mn in gravitational interaction, and let O be their centre of mass. For every body Pi, we denote ui = OPi. In the barycentric reference frame with origin O, Newton’s equations of motion form a differential system of order 6(n + 1) and can be written in Hamiltonian form using the canonical coordinates with Hamiltonian (1)where Δij = ∥ uiuj, and is the constant of gravitation. The reduction of the centre of mass is achieved by using the canonical heliocentric variables of Poincaré (Poincaré 1905; Laskar & Robutel 1995), defined as This Hamiltonian can then be split into an integrable part, H0, and a perturbation, H1, (4)with (5)and (6)The integrable part, H0, is the Hamiltonian of a sum of disjoint Kepler problems of a single planet of mass mi around a fixed Sun of mass m0. A set of adapted variables for H0 will thus be given by the elliptical elements, (ak,ek,ik,λk,ϖkk), where ak is the semi-major axis, ek the eccentricity, ik the inclination, λk the mean longitude, ϖk the longitude of the perihelion, and Ωk the longitude of the node. They are defined as the elliptical elements associated to the Hamiltonian (7)with .

2.2. Angular momentum deficit (AMD)

The total linear momentum of the system is null in the barycentric reference frame (8)Let G be the total angular momentum. Its expression is not changed by the linear symplectic change of variable . We have thus (9)When expressed in heliocentric variables, the angular momentum is thus the sum of the angular momentum of the Keplerian problems of the unperturbed Hamiltonian H0. In particular, if the angular momentum direction is assumed to be the axis z, the norm of the angular momentum is (10)where . For such a system, the AMD is defined as the difference between the norm of the angular momentum of a coplanar and circular system with the same semi-major axis values and the norm of the angular momentum (G), i.e. (Laskar 1997, 2000) (11)

2.3. AMD and collision of orbits

The instabilities of a planetary system often result in a modification of its architecture. A planet can be ejected from the system or can fall into the star; in both cases this results in a loss of AMD for the system. The outcome of the AMD after a planetary collisions is less trivial and needs to be computed. Assume that among our n + 1 bodies, the (totally inelastic) collision of two bodies of masses m1 and m2, and orbits occurs, forming a new body . During this collision we assume that the other bodies are not affected. The mass is conserved (12)and the linear momentum in the barycentric reference frame is conserved so , and also (13)On the other hand, r3 = r1 = r2 at the time of the collision, so the angular momentum is also conserved (14)It should be noted that the transformation of the orbits during the collision is perfectly defined by Eqs. (12), (13). The problem which remains is to compute the evolution of the elliptical elements during the collision.

2.3.1. Energy evolution during collision

Just before the collision, we assume that the orbits and are elliptical heliocentric orbits. At the time of the collision, only these two bodies are involved and the other bodies are not affected. The evolution of the orbits are thus given by the conservation laws Eqs. (12), (13). The Keplerian energy of each particle is (15)At collision, we have r1 = r2 = r3 = r; we have thus the conservation of the potential energy (16)The change of Keplerian energy is thus given by the change of kinetic energy (17)that is, with Eqs. (12), (13), (18)Thus, the Keplerian energy of the system decreases during collision. Part of the kinetic energy is dispelled during the collision. As expected, there is no loss of energy when , that is, as m1m2 ≠ 0, when . As an immediate consequence of the decrease of energy during the collision, we have (19)with (20)

2.3.2. AMD evolution during collision

Let . As f′(x) < 0 and f′′(x) > 0, we have, as f is decreasing and convex, (21)thus (22)During the collision, the angular momentum is conserved Eq. (14), and so is the conservation of its normal component, that is (23)We deduce that in all circumstances we have a decrease of the angular momentum deficit during the collision, that is (24)with (25)The equality can hold in Eq. (24) only if m1 = 0, m2 = 0, or a1 = a2 and , that is when one of the bodies is massless, or when the two bodies are on the same orbit and at the same position (at the time of the collision, we also have r1 = r2).

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

3. AMD-stability

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

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

3.1. Collisional condition

We assume that a,a,m,m are non-zero. Denoting α = a/a, ≳ = m/m, the system in Eq. (27) becomes (28)(29)with , and where C(e,e′) = C/ Λ′ is called the relative AMD. As e and e are planetary eccentricities, we also have (30)

thumbnail Fig. 1

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

Open with DEXTER

The set of collisional conditions (Eqs. (28)–(30)) can be solved using Lagrange multipliers. We are looking for the minimum value of the relative AMD C(e,e′) (29) for which the collision condition (28) is satisfied. These conditions are visualised in the (e,e′) plane in Fig. 1. The collision condition (28) corresponds to the segment AB of Fig. 1. The domain of collisions is limited by the conditions (30). For e = 0, we have , and the intersection of the collision line with the axis e′ = 0 is obtained for e0 = 1 /α − 1. This value can be greater or smaller than 1 depending on the value of α. We have thus the different cases (Fig. 1) and the limit case α = 1 / 2, for which e0 = 1. In all cases, the Lagrange multipliers condition is written (33)which gives (34)This relation allows e to be eliminated in the collision condition (28), which becomes an equation in the single variable e, and parameters (α, ≳): (35)Here F(e,a, ≳) is properly defined for (e,a, ≳) in the domain De,α, defined by e ∈ [ 0,1 ] ∈ ] 0,1 ] , ≳ ∈ ] 0, + ∞ [, as in this domain, 1 − e2 + ≳ 2e2/α> 0. We also have (36)Thus, on the domain De,α, and F(e,α, ≳) is strictly increasing with respect to e for e ∈ [ 0,1 ]. Moreover, as 0 <α< 1, (37)and (38)The equation of collision (35) thus always has a single solution ec in the interval ] 0,min(1,e0) [. This ensures that this critical value of e will fulfil the condition (30). The corresponding value of the relative AMD is then obtained through (29).

3.2. Critical AMD C

We have thus demonstrated that for a given pair of ratios of semi-major axes and masses, (α, ≳), there is always a unique critical value Cc(α, ≳) of the relative AMD C= C/ Λ′ which defines the AMD-stability. The system of two planets is AMD-stable if and only if (39)The value of the critical AMD Cc(α, ≳) is obtained by computing first the critical eccentricity ec(α, ≳) which is the unique solution of the collisional Eq. (35) in the interval [ 0,1 ]. The critical AMD is then (Eq. (29)), where the critical value is obtained from ec through Eq. (28). It is important to note that the critical AMD, and thus the AMD-stability condition, depends only on (α, ≳).

Table 1

Special values of Cc(α, ≳).

3.3. Behaviour of the critical AMD

We now analyse the general properties of the critical AMD function Cc(α, ≳). As , we can apply the implicit function theorem to the domain De,α,, which then ensures that the solution of the collision Eq. (35), ec(α, ≳), is a continuous function of (and even analytic for ≳ ∈ ] 0, + ∞ [). Moreover, on De,α,,(40)We also have(41)and ec(α, ≳) is a decreasing function of . For any given values of the semi-major axes ratio α, and masses , we can thus find the critical value Cc(α, ≳) which allows for a collision (39). For the critical value Cc(α, ≳), a single solution corresponds to the tangency condition (Fig. 1), and this solution is obtained at the critical value ec(α, ≳) for the eccentricity of the orbit . The values of the critical relative AMD Cc(α, ≳) are plotted in Fig. 2 versus α, for different values of . Deriving Eq. (29)with respect to , one obtains (42)Using the two relations (Eqs. (28), (34)), this reduces to (43)Thus Cc(α, ≳) strictly increases with . In the same way (44)and Cc(α, ≳) decreases with α.

thumbnail Fig. 2

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

Open with DEXTER

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

A development of Cc(α, ≳) for α − → 1 can also be obtained (see Appendix C). With η = 1 − α, we have (45)We now present two applications of the AMD-stability, the formation of planetary systems, and a classification of observed planetary systems.

4. AMD model of planet formation

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

4.1. Collisional evolution

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

In the accretion process, we consider a planetesimal of semi-major axis a and its immediate neighbour, defined as the planetesimal with semi-major axis a, such that there is no other planetesimal with semi-major axis between a and a. In this case we can assume that α is close to 1 and, as explained in Appendix C, we use an approximation of the critical AMD value Cc(α,γ): (46)where δa = a′ − a and k( ≳ ) = ≳ / (2( ≳ + 1)).

4.2. AMD-stable planetary distribution

In this section, we search for the laws followed by the planetary distribution of a model formed following the above assumptions. We thus start from an arbitrary distribution of mass of planetesimals ρ(a), and let the system evolve under the previous rules. We search for the condition of AMD-stable planetary systems, obtained by random accretion of planetesimals. This condition requires that the final AMD value cannot allow for planetary crossing among the planets. Let C be the value of the AMD at the end of the accretion process. If we consider the planetesimal of semi-major axis a, its mass will continue to increase by accretion with a body of semi-major axis a′>a, as long as (47)The initial linear density of mass is ρ(a). As a is the closest neighbour to a, we can assume that all the planetesimals initially between a and a have been absorbed by the two bodies of mass m(a) and m(a′). At first order with respect to δa/a, we have thus (48)and from (47) at the limiting case we have (49)where . We have thus (50)and from Eq. (48) (51)

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

Using the previous relations, we can compute the resulting systems for various initial mass distribution, in particular for ρ(a) = ζap. From Eq. (51), we obtain for two consecutive planets (52)and from Eq. (46), as limα − → 1 ≳ = 1, we thus have in all cases. Relation (50) can be rewritten (53)where δn = 1 is the increment from planet a to a. By integration, this difference relation becomes for , (54)For , we obtain (55)In particular, for p = 0 (constant distribution and ρ(a) = ζ), from Eq. (53), we have (56)For the masses, from Eq. (51), we have for large n (or if a0 is small) (57)For p = − 3 / 2, we find a power law similar to the Titius-Bode law for planetary distribution2. The different expressions deduced from this model of planetary accretion are listed in Table 2.

Table 2

Planetary distribution corresponding to different initial mass distribution.

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

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

5. AMD-classification of planetary systems

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

5.1. Sample studied and methods of computation

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

thumbnail Fig. 3

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

Open with DEXTER

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

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

Table 3

Result of the analysis split in function of the multiplicity of the system.

thumbnail Fig. 4

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

Open with DEXTER

5.2. Propagation of uncertainties

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

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

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

thumbnail Fig. 5

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

Open with DEXTER

thumbnail Fig. 6

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

Open with DEXTER

5.3. AMD-stable systems

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

5.3.1. Strong AMD-stable systems

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

5.3.2. Weak AMD-stable systems

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

5.4. AMD-unstable systems

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

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

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

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

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

5.4.2. Resonant systems

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

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

5.4.3. Concentration around MMR

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

thumbnail Fig. 7

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

Open with DEXTER

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

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

Let us denote the probability density of the catalogue period ratios f and the associated cumulative distribution . Let us consider an interval Ix = (x,x + Δx), if we assume 0, the probability for a ratio r to be in Ix is (64)Therefore, the probability that more than k0 pairs out of N pairs drawn randomly from the distribution are in Ix is (65)This is the probability of a binomial trial. From now on, N designates the number of period ratios in Ru.

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

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

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

thumbnail Fig. 8

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

Open with DEXTER

thumbnail Fig. 9

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

Open with DEXTER

thumbnail Fig. 10

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

Open with DEXTER

thumbnail Fig. 11

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

Open with DEXTER

thumbnail Fig. 12

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

Open with DEXTER

6. Conclusion

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

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

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

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

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


1

Here practical stability means stability over a very long time compared to the expected life of the central star.

2

The value p = − 3 / 2 corresponds to a surface density proportional to r− 5 / 2, which is different from the − 3 / 2 surface density exponent of the minimum mass solar nebula (Weidenschilling 1977). For the solar system, (Laskar 2000) found p = 0 to be the best fitting value; it corresponds to a surface density proportional to r-1.

3

This is equivalent to require that all pairs are AMD-stable.

References

Appendix A: Intersection of planar orbits

thumbnail Fig. A.1

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

Open with DEXTER

In this section, we present an efficient algorithm for the computation of the intersection of two elliptical orbits in the plane, following (Albouy 2002). Let us consider an elliptical orbit defined by (μ,r,) and let G = r be the angular momentum per unit of mass. The Laplace vector (A.1)is an integral of the motion with coordinates (ecosω,esinω). One has (A.2)where p = a(1 − e2) is the parameter of the ellipse. Let r = (x,y) in the plane. We can consider the ellipse in three-dimensional space (x,y,r) as the intersection of the cone (A.3)with the plane defined by Eq. (A.2), that is (A.4)If we consider now two orbits and . Their intersection is easily obtained as the intersection of the line of intersection of the two planes (A.5)with the cone of equation r2 = x2 + y2. Depending on the initial conditions, if and are distinct, we will get either 0, 1, or 2 solutions.

Appendix B: AMD in the averaged equations

In this section we show that the AMD conserves the same form in the averaged planetary Hamiltonian at all orders. More generally, this is true for any integral of H which does not depend on the longitude λk. Let (B.1)be a perturbed Hamiltonian system, and let K,J,φ) be a first integral of H,λ,J,φ) (such that { K,H } = 0, where { ·,· } is the usual Poisson bracket), and independent of λ. In addition, let (B.2)be a formal averaging of H with respect to λ. If S,λ,J,φ) is a generator defined as below (Eqs. (B.4), (B.6), (B.7)), such that H is independent of λ. Then, K is an integral of H, i.e. (B.3)The generator S = εS1 + ε2S2 + ε3S3 +··· is obtained formally through the identification order by order (B.4)For any function G,λ,J,φ), let (B.5)be the average of G over all the angles λk. For each n ≥ 1, the Sn is obtained through the resolution of an equation of the form (B.6)where n belongs to ℒ(H0,H1,S1,...,Sn − 1), the Lie algebra generated by (H0,H1,S1,...,Sn − 1). will be the averaged part of n, ⟨ ℛnλ, and Sn is obtained by solving the homological equation (B.7)We show by recurrence that { K,Sn } = 0 for all n ≥ 1. First, we note that as { K,H0 } = 0, we also have { K,H1 } = 0. As K is independent of λk, we also have for all G(B.8)This can be seen by formal expansion of G in a Fourier series G = ∑ gkeik,λ. We have thus Gλ = g0. Let us assume now that { K,Sk } = 0 for all kn. As n + 1 ∈ ℒ(H0,H1,S1,...,Sn), we also have { K,n + 1 } = 0, and from Eq. (B.8) { K, ⟨ ℛn + 1λ } = 0 and thus (B.9)Using the Jacobi identity, as { F,H0 } = 0, we have (B.10)The solution of the homological Eq. (B.7) is unique up to a term independent of λ. But as ⟨ { K,Sn + 1 } ⟩ λ = 0, then the only possible solution for Eq. (B.10) is (B.11)In the same way, as , we have . Thus { K, { H0,S1 } } = 0, by the Jacobi identity, { H0, { K,S1 } } = 0, and as previously, { K,S1 } = 0. Our recurrence is thus complete and it follows immediately that { K,H′ } = 0.

Appendix C: Special values of C

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

Appendix C.1: Asymptotic value of C for ≳ → 0

We have shown that ec(α, ≳) is a differentiable function of , which is monotonic (41) and bounded (ec(a, ≳ ) ∈ [ 0,1 ]). The limit exists for all α ∈ ] 0,1 ]. If ec(α, ≳) is a solution of equation(35), it will also be a solution of the following cubic equation (in e), which is directly obtained from (35) by squaring, multiplication, and simplification by α(1 + e):(C.1)As ec(α, ≳), is a continuous function of , when ≳ → 0 the limits e0c(α) will satisfy the limit equation (C.2)with solutions e0 = 1 /α − 1 and e1 = 1. Depending on α, several cases are treated:

α< 1 / 2:

we have then e0> 1, and the only possibility for ec(α,0) is e1 = 1; as it is the only root of Eq. (C.2) which belongs to [ 0,1 ] , we have thus (C.3)We have then (C.4)In order to study the behaviour of ec(α, ≳) in the vicinity of ≳ = 0, we can differentiate K(ec(α, ≳ ), ≳ ) = 0 twice, which gives (C.5)thus ec(α, ≳) decreases with respect to at ≳ = 0.

The second order development of Cc gives (C.6)

α> 1 / 2:

in this case, e0< 1. As ec(α, ≳ ) ∈ ] 0,e0 [, we have ec(α,0) ∈ [ 0,e0 ], which gives as the unique possibility (C.7)with (C.8)By setting ≳ = 0 in Eq. (41), we also obtain (C.9)The development of Cc gives (C.10)

α = 1 / 2:

in this case, e0 = 1, and the only solution is (C.11)and (C.12)Moreover, equation (C.1) becomes 22e2(3 − e) = (1 − e)3. We obtain thus the asymptotic value for ec(1 / 2, ≳) when ≳ → 0 as (C.13)For α = 1 / 2, the development of Cc in γ contains non-polynomial terms in giving (C.14)

Appendix C.2: Asymptotic value of C for ≳ → + ∞

This case is more simple. If ec(α, ≳) is a solution of Eq. (35), then it will also be a solution of Eq. (C.1), and thus of (C.15)As ec(α, ≳) is monotonic and bounded, it has a limit when ≳ → + ∞, which will verify the limit equation (C.15), when ≳ → + ∞, that is (C.16)As 0 <α< 1, the only solution is e = 0, and thus (C.17)and (C.18)and more precisely (C.19)

Appendix C.3: Study of C for ≳ = 1 and

For ≳ = 1 or , we can also obtain simple expressions for Cc(α, ≳). Indeed, If ≳ = 1, K(e,α, ≳) factorises in (1 − α)(1 + e)(αe2 − 2e + 1 − α) and we have a single solution for ec in the interval [ 0,1 ], (C.20)and (C.21)For , the cubic Eq. (C.1) reduces to (C.22)with the single solution (C.23)and (C.24)

Appendix C.4: C for α − → 1

Let us denote η = 1 − α. Equation (C.1)can be developed in η, (C.25)The zeroth and first orders of Eq. (C.25)imply that ec must go to zero; moreover, it scales with η. We write ec(η, ≳ ) = κ( ≳ )η + o(η). We inject this expression in Eq. (C.25)and keep the second order in η(C.26)We keep the solution that is positive and continuous in and we have (C.27)If we now compute Cc developed for α − → 1 we have (C.28)

All Tables

Table 1

Special values of Cc(α, ≳).

Table 2

Planetary distribution corresponding to different initial mass distribution.

Table 3

Result of the analysis split in function of the multiplicity of the system.

All Figures

thumbnail Fig. 1

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

Open with DEXTER
In the text
thumbnail Fig. 2

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

Open with DEXTER
In the text
thumbnail Fig. 3

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

Open with DEXTER
In the text
thumbnail Fig. 4

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

Open with DEXTER
In the text
thumbnail Fig. 5

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

Open with DEXTER
In the text
thumbnail Fig. 6

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

Open with DEXTER
In the text
thumbnail Fig. 7

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

Open with DEXTER
In the text
thumbnail Fig. 8

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

Open with DEXTER
In the text
thumbnail Fig. 9

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

Open with DEXTER
In the text
thumbnail Fig. 10

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

Open with DEXTER
In the text
thumbnail Fig. 11

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

Open with DEXTER
In the text
thumbnail Fig. 12

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

Open with DEXTER
In the text
thumbnail Fig. A.1

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

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.